Actual source code: ex37.c

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

  3: #include <petscmat.h>

  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;

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

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

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

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

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

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

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

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

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

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

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;

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

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

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

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

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

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

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

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

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

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

169: /*TEST

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}

175: TEST*/