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