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 PetscCall(PetscInitialize(&argc,&args,(char*)0,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.; v[1] = 2.; v[2] = -1.; 40 for (i=1; i<n-1; i++) { 41 midx[2] = i-1; midx[1] = i; midx[0] = i+1; 42 PetscCall(MatSetValues(C,1,&i,3,midx,v,INSERT_VALUES)); 43 } 44 i = 0; midx[0] = 0; midx[1] = 1; 45 v[0] = 2.0; v[1] = -1.; 46 PetscCall(MatSetValues(C,1,&i,2,midx,v,INSERT_VALUES)); 47 i = n-1; midx[0] = n-2; midx[1] = n-1; 48 v[0] = -1.0; v[1] = 2.; 49 PetscCall(MatSetValues(C,1,&i,2,midx,v,INSERT_VALUES)); 50 51 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 52 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 53 PetscCall(MatView(C,NULL)); 54 PetscCall(MatViewFromOptions(C,NULL,"-view")); 55 56 /* test matduplicate */ 57 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&A)); 58 PetscCall(PetscObjectSetName((PetscObject)A,"duplicate_copy")); 59 PetscCall(MatViewFromOptions(A,NULL,"-view")); 60 PetscCall(MatEqual(A,C,&flg)); 61 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDuplicate(C,MAT_COPY_VALUES,): Matrices are NOT equal"); 62 PetscCall(MatDestroy(&A)); 63 64 /* test matrices with different nonzero patterns - Note: A is created with different nonzero pattern of C! */ 65 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 66 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n)); 67 PetscCall(MatSetFromOptions(A)); 68 PetscCall(MatSetUp(A)); 69 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 70 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 71 72 PetscCall(MatCopy(C,A,DIFFERENT_NONZERO_PATTERN)); 73 PetscCall(PetscObjectSetName((PetscObject)A,"copy_diffnnz")); 74 PetscCall(MatViewFromOptions(A,NULL,"-view")); 75 PetscCall(MatEqual(A,C,&flg)); 76 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatCopy(C,A,DIFFERENT_NONZERO_PATTERN): Matrices are NOT equal"); 77 78 /* test matrices with same nonzero pattern */ 79 PetscCall(MatDestroy(&A)); 80 PetscCall(MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&A)); 81 PetscCall(MatCopy(C,A,SAME_NONZERO_PATTERN)); 82 PetscCall(PetscObjectSetName((PetscObject)A,"copy_samennz")); 83 PetscCall(MatViewFromOptions(A,NULL,"-view")); 84 PetscCall(MatEqual(A,C,&flg)); 85 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatCopy(C,A,SAME_NONZERO_PATTERN): Matrices are NOT equal"); 86 87 /* test subset nonzero pattern */ 88 PetscCall(MatCopy(C,A,SUBSET_NONZERO_PATTERN)); 89 PetscCall(PetscObjectSetName((PetscObject)A,"copy_subnnz")); 90 PetscCall(MatViewFromOptions(A,NULL,"-view")); 91 PetscCall(MatEqual(A,C,&flg)); 92 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatCopy(C,A,SUBSET_NONZERO_PATTERN): Matrices are NOT equal"); 93 94 /* Test MatCopy on a matrix obtained after MatConvert from AIJ 95 see https://lists.mcs.anl.gov/pipermail/petsc-dev/2019-April/024289.html */ 96 PetscCall(MatHasCongruentLayouts(C,&flg)); 97 if (flg) { 98 Mat Cs,Cse; 99 MatType Ctype,Cstype; 100 101 PetscCall(MatGetType(C,&Ctype)); 102 PetscCall(MatTranspose(C,MAT_INITIAL_MATRIX,&Cs)); 103 PetscCall(MatAXPY(Cs,1.0,C,DIFFERENT_NONZERO_PATTERN)); 104 PetscCall(MatConvert(Cs,MATAIJ,MAT_INPLACE_MATRIX,&Cs)); 105 PetscCall(MatSetOption(Cs,MAT_SYMMETRIC,PETSC_TRUE)); 106 PetscCall(MatGetType(Cs,&Cstype)); 107 108 PetscCall(PetscObjectSetName((PetscObject)Cs,"symm_initial")); 109 PetscCall(MatViewFromOptions(Cs,NULL,"-view")); 110 111 PetscCall(MatConvert(Cs,Ctype,MAT_INITIAL_MATRIX,&Cse)); 112 PetscCall(PetscObjectSetName((PetscObject)Cse,"symm_conv_init")); 113 PetscCall(MatViewFromOptions(Cse,NULL,"-view")); 114 PetscCall(MatMultEqual(Cs,Cse,5,&flg)); 115 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatConvert MAT_INITIAL_MATRIX %s -> %s: Matrices are NOT multequal",Ctype,Cstype); 116 117 PetscCall(MatConvert(Cs,Ctype,MAT_REUSE_MATRIX,&Cse)); 118 PetscCall(PetscObjectSetName((PetscObject)Cse,"symm_conv_reuse")); 119 PetscCall(MatViewFromOptions(Cse,NULL,"-view")); 120 PetscCall(MatMultEqual(Cs,Cse,5,&flg)); 121 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatConvert MAT_REUSE_MATRIX %s -> %s: Matrices are NOT multequal",Ctype,Cstype); 122 123 PetscCall(MatCopy(Cs,Cse,SAME_NONZERO_PATTERN)); 124 PetscCall(PetscObjectSetName((PetscObject)Cse,"symm_conv_copy_samennz")); 125 PetscCall(MatViewFromOptions(Cse,NULL,"-view")); 126 PetscCall(MatMultEqual(Cs,Cse,5,&flg)); 127 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatCopy(...SAME_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal",Ctype,Cstype); 128 129 PetscCall(MatCopy(Cs,Cse,SUBSET_NONZERO_PATTERN)); 130 PetscCall(PetscObjectSetName((PetscObject)Cse,"symm_conv_copy_subnnz")); 131 PetscCall(MatViewFromOptions(Cse,NULL,"-view")); 132 PetscCall(MatMultEqual(Cs,Cse,5,&flg)); 133 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatCopy(...SUBSET_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal",Ctype,Cstype); 134 135 PetscCall(MatCopy(Cs,Cse,DIFFERENT_NONZERO_PATTERN)); 136 PetscCall(PetscObjectSetName((PetscObject)Cse,"symm_conv_copy_diffnnz")); 137 PetscCall(MatViewFromOptions(Cse,NULL,"-view")); 138 PetscCall(MatMultEqual(Cs,Cse,5,&flg)); 139 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatCopy(...DIFFERENT_NONZERO_PATTERN) %s -> %s: Matrices are NOT multequal",Ctype,Cstype); 140 141 PetscCall(MatDestroy(&Cse)); 142 PetscCall(MatDestroy(&Cs)); 143 } 144 145 /* test MatStore/RetrieveValues() */ 146 if (isAIJ) { 147 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE)); 148 PetscCall(MatStoreValues(A)); 149 PetscCall(MatZeroEntries(A)); 150 PetscCall(MatRetrieveValues(A)); 151 } 152 153 PetscCall(MatDestroy(&C)); 154 PetscCall(MatDestroy(&A)); 155 PetscCall(PetscFinalize()); 156 return 0; 157 } 158 159 /*TEST 160 161 testset: 162 nsize: {{1 2}separate output} 163 args: -view ::ascii_info -mat_type {{aij baij sbaij mpiaij mpibaij mpisbaij}separate output} -mat_block_size {{1 2}separate output} 164 165 TEST*/ 166