1 static char help[] = "Tests MatCopy() and MatStore/RetrieveValues().\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **args) 6 { 7 Mat C, A; 8 PetscInt i, n = 10, midx[3], bs = 1; 9 PetscScalar v[3]; 10 PetscBool flg, isAIJ; 11 MatType type; 12 PetscMPIInt size; 13 14 PetscFunctionBeginUser; 15 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 16 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 17 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL)); 19 20 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 21 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, n, n)); 22 PetscCall(MatSetType(C, MATAIJ)); 23 PetscCall(MatSetFromOptions(C)); 24 PetscCall(PetscObjectSetName((PetscObject)C, "initial")); 25 26 PetscCall(MatGetType(C, &type)); 27 if (size == 1) { 28 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATSEQAIJ, &isAIJ)); 29 } else { 30 PetscCall(PetscObjectTypeCompare((PetscObject)C, MATMPIAIJ, &isAIJ)); 31 } 32 PetscCall(MatSeqAIJSetPreallocation(C, 3, NULL)); 33 PetscCall(MatMPIAIJSetPreallocation(C, 3, NULL, 3, NULL)); 34 PetscCall(MatSeqBAIJSetPreallocation(C, bs, 3, NULL)); 35 PetscCall(MatMPIBAIJSetPreallocation(C, bs, 3, NULL, 3, NULL)); 36 PetscCall(MatSeqSBAIJSetPreallocation(C, bs, 3, NULL)); 37 PetscCall(MatMPISBAIJSetPreallocation(C, bs, 3, NULL, 3, NULL)); 38 39 v[0] = -1.; 40 v[1] = 2.; 41 v[2] = -1.; 42 for (i = 1; i < n - 1; i++) { 43 midx[2] = i - 1; 44 midx[1] = i; 45 midx[0] = i + 1; 46 PetscCall(MatSetValues(C, 1, &i, 3, midx, v, INSERT_VALUES)); 47 } 48 i = 0; 49 midx[0] = 0; 50 midx[1] = 1; 51 v[0] = 2.0; 52 v[1] = -1.; 53 PetscCall(MatSetValues(C, 1, &i, 2, midx, v, INSERT_VALUES)); 54 i = n - 1; 55 midx[0] = n - 2; 56 midx[1] = n - 1; 57 v[0] = -1.0; 58 v[1] = 2.; 59 PetscCall(MatSetValues(C, 1, &i, 2, midx, v, INSERT_VALUES)); 60 61 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 62 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 63 PetscCall(MatView(C, NULL)); 64 PetscCall(MatViewFromOptions(C, NULL, "-view")); 65 66 /* test matduplicate */ 67 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &A)); 68 PetscCall(PetscObjectSetName((PetscObject)A, "duplicate_copy")); 69 PetscCall(MatViewFromOptions(A, NULL, "-view")); 70 PetscCall(MatEqual(A, C, &flg)); 71 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatDuplicate(C,MAT_COPY_VALUES,): Matrices are NOT equal"); 72 PetscCall(MatDestroy(&A)); 73 74 /* test matrices with different nonzero patterns - Note: A is created with different nonzero pattern of C! */ 75 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 76 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)); 77 PetscCall(MatSetFromOptions(A)); 78 PetscCall(MatSetUp(A)); 79 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 80 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 81 82 PetscCall(MatCopy(C, A, DIFFERENT_NONZERO_PATTERN)); 83 PetscCall(PetscObjectSetName((PetscObject)A, "copy_diffnnz")); 84 PetscCall(MatViewFromOptions(A, NULL, "-view")); 85 PetscCall(MatEqual(A, C, &flg)); 86 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatCopy(C,A,DIFFERENT_NONZERO_PATTERN): Matrices are NOT equal"); 87 88 /* test matrices with same nonzero pattern */ 89 PetscCall(MatDestroy(&A)); 90 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &A)); 91 PetscCall(MatCopy(C, A, SAME_NONZERO_PATTERN)); 92 PetscCall(PetscObjectSetName((PetscObject)A, "copy_samennz")); 93 PetscCall(MatViewFromOptions(A, NULL, "-view")); 94 PetscCall(MatEqual(A, C, &flg)); 95 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatCopy(C,A,SAME_NONZERO_PATTERN): Matrices are NOT equal"); 96 97 /* test subset nonzero pattern */ 98 PetscCall(MatCopy(C, A, SUBSET_NONZERO_PATTERN)); 99 PetscCall(PetscObjectSetName((PetscObject)A, "copy_subnnz")); 100 PetscCall(MatViewFromOptions(A, NULL, "-view")); 101 PetscCall(MatEqual(A, C, &flg)); 102 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatCopy(C,A,SUBSET_NONZERO_PATTERN): Matrices are NOT equal"); 103 104 /* Test MatCopy on a matrix obtained after MatConvert from AIJ 105 see https://lists.mcs.anl.gov/pipermail/petsc-dev/2019-April/024289.html */ 106 PetscCall(MatHasCongruentLayouts(C, &flg)); 107 if (flg) { 108 Mat Cs, Cse; 109 MatType Ctype, Cstype; 110 111 PetscCall(MatGetType(C, &Ctype)); 112 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Cs)); 113 PetscCall(MatAXPY(Cs, 1.0, C, DIFFERENT_NONZERO_PATTERN)); 114 PetscCall(MatConvert(Cs, MATAIJ, MAT_INPLACE_MATRIX, &Cs)); 115 PetscCall(MatSetOption(Cs, MAT_SYMMETRIC, PETSC_TRUE)); 116 PetscCall(MatGetType(Cs, &Cstype)); 117 118 PetscCall(PetscObjectSetName((PetscObject)Cs, "symm_initial")); 119 PetscCall(MatViewFromOptions(Cs, NULL, "-view")); 120 121 PetscCall(MatConvert(Cs, Ctype, MAT_INITIAL_MATRIX, &Cse)); 122 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_init")); 123 PetscCall(MatViewFromOptions(Cse, NULL, "-view")); 124 PetscCall(MatMultEqual(Cs, Cse, 5, &flg)); 125 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatConvert MAT_INITIAL_MATRIX %s -> %s: Matrices are NOT multequal", Ctype, Cstype); 126 127 PetscCall(MatConvert(Cs, Ctype, MAT_REUSE_MATRIX, &Cse)); 128 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_reuse")); 129 PetscCall(MatViewFromOptions(Cse, NULL, "-view")); 130 PetscCall(MatMultEqual(Cs, Cse, 5, &flg)); 131 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatConvert MAT_REUSE_MATRIX %s -> %s: Matrices are NOT multequal", Ctype, Cstype); 132 133 PetscCall(MatCopy(Cs, Cse, SAME_NONZERO_PATTERN)); 134 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_copy_samennz")); 135 PetscCall(MatViewFromOptions(Cse, NULL, "-view")); 136 PetscCall(MatMultEqual(Cs, Cse, 5, &flg)); 137 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatCopy(...SAME_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal", Ctype, Cstype); 138 139 PetscCall(MatCopy(Cs, Cse, SUBSET_NONZERO_PATTERN)); 140 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_copy_subnnz")); 141 PetscCall(MatViewFromOptions(Cse, NULL, "-view")); 142 PetscCall(MatMultEqual(Cs, Cse, 5, &flg)); 143 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatCopy(...SUBSET_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal", Ctype, Cstype); 144 145 PetscCall(MatCopy(Cs, Cse, DIFFERENT_NONZERO_PATTERN)); 146 PetscCall(PetscObjectSetName((PetscObject)Cse, "symm_conv_copy_diffnnz")); 147 PetscCall(MatViewFromOptions(Cse, NULL, "-view")); 148 PetscCall(MatMultEqual(Cs, Cse, 5, &flg)); 149 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatCopy(...DIFFERENT_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal", Ctype, Cstype); 150 151 PetscCall(MatDestroy(&Cse)); 152 PetscCall(MatDestroy(&Cs)); 153 } 154 155 /* test MatStore/RetrieveValues() */ 156 if (isAIJ) { 157 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE)); 158 PetscCall(MatStoreValues(A)); 159 PetscCall(MatZeroEntries(A)); 160 PetscCall(MatRetrieveValues(A)); 161 } 162 163 PetscCall(MatDestroy(&C)); 164 PetscCall(MatDestroy(&A)); 165 PetscCall(PetscFinalize()); 166 return 0; 167 } 168 169 /*TEST 170 171 testset: 172 nsize: {{1 2}separate output} 173 args: -view ::ascii_info -mat_type {{aij baij sbaij mpiaij mpibaij mpisbaij}separate output} -mat_block_size {{1 2}separate output} 174 175 TEST*/ 176