Actual source code: ex37.c


  2: static char help[] = "Tests MatCopy() and MatStore/RetrieveValues().\n\n";

  4: #include <petscmat.h>

  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;

 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));

 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"));

 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));

 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));

 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"));

 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));

 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));

 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");

 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");

 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");

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;

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));

119:     PetscCall(PetscObjectSetName((PetscObject)Cs, "symm_initial"));
120:     PetscCall(MatViewFromOptions(Cs, NULL, "-view"));

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);

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);

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);

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);

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);

152:     PetscCall(MatDestroy(&Cse));
153:     PetscCall(MatDestroy(&Cs));
154:   }

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:   }

164:   PetscCall(MatDestroy(&C));
165:   PetscCall(MatDestroy(&A));
166:   PetscCall(PetscFinalize());
167:   return 0;
168: }

170: /*TEST

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}

176: TEST*/