Actual source code: ex54.c

  1: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel AIJ and BAIJ formats.\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat          E, A, B, Bt, *submatA, *submatB;
  8:   PetscInt     bs = 1, m = 11, ov = 1, i, j, k, *rows, *cols, nd = 5, *idx, rstart, rend, sz, mm, nn, M, N, Mbs;
  9:   PetscMPIInt  size, rank;
 10:   PetscScalar *vals, rval;
 11:   IS          *is1, *is2;
 12:   PetscRandom  rdm;
 13:   Vec          xx, s1, s2;
 14:   PetscReal    s1norm, s2norm, rnorm, tol = 100 * PETSC_SMALL;
 15:   PetscBool    flg, test_nd0 = PETSC_FALSE, emptynd;

 17:   PetscFunctionBeginUser;
 18:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 19:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 20:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 22:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
 23:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL));
 24:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
 25:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
 26:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nd0", &test_nd0, NULL));

 28:   /* Create a AIJ matrix A */
 29:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 30:   PetscCall(MatSetSizes(A, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE));
 31:   PetscCall(MatSetType(A, MATAIJ));
 32:   PetscCall(MatSetBlockSize(A, bs));
 33:   PetscCall(MatSeqAIJSetPreallocation(A, PETSC_DEFAULT, NULL));
 34:   PetscCall(MatMPIAIJSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
 35:   PetscCall(MatSetFromOptions(A));
 36:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));

 38:   /* Create a BAIJ matrix B */
 39:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
 40:   PetscCall(MatSetSizes(B, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE));
 41:   PetscCall(MatSetType(B, MATBAIJ));
 42:   PetscCall(MatSeqBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL));
 43:   PetscCall(MatMPIBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
 44:   PetscCall(MatSetFromOptions(B));
 45:   PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));

 47:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
 48:   PetscCall(PetscRandomSetFromOptions(rdm));

 50:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 51:   PetscCall(MatGetSize(A, &M, &N));
 52:   Mbs = M / bs;

 54:   PetscCall(PetscMalloc1(bs, &rows));
 55:   PetscCall(PetscMalloc1(bs, &cols));
 56:   PetscCall(PetscMalloc1(bs * bs, &vals));
 57:   PetscCall(PetscMalloc1(M, &idx));

 59:   /* Now set blocks of values */
 60:   for (i = 0; i < 40 * bs; i++) {
 61:     PetscInt nr = 1, nc = 1;
 62:     PetscCall(PetscRandomGetValue(rdm, &rval));
 63:     cols[0] = bs * (int)(PetscRealPart(rval) * Mbs);
 64:     PetscCall(PetscRandomGetValue(rdm, &rval));
 65:     rows[0] = rstart + bs * (int)(PetscRealPart(rval) * m);
 66:     for (j = 1; j < bs; j++) {
 67:       PetscCall(PetscRandomGetValue(rdm, &rval));
 68:       if (PetscRealPart(rval) > .5) rows[nr++] = rows[0] + j - 1;
 69:     }
 70:     for (j = 1; j < bs; j++) {
 71:       PetscCall(PetscRandomGetValue(rdm, &rval));
 72:       if (PetscRealPart(rval) > .5) cols[nc++] = cols[0] + j - 1;
 73:     }

 75:     for (j = 0; j < nr * nc; j++) {
 76:       PetscCall(PetscRandomGetValue(rdm, &rval));
 77:       vals[j] = rval;
 78:     }
 79:     PetscCall(MatSetValues(A, nr, rows, nc, cols, vals, ADD_VALUES));
 80:     PetscCall(MatSetValues(B, nr, rows, nc, cols, vals, ADD_VALUES));
 81:   }
 82:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 83:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 84:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 85:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

 87:   /* Test MatConvert_MPIAIJ_MPI(S)BAIJ handles incompletely filled blocks */
 88:   PetscCall(MatConvert(A, MATBAIJ, MAT_INITIAL_MATRIX, &E));
 89:   PetscCall(MatDestroy(&E));
 90:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Bt));
 91:   PetscCall(MatAXPY(Bt, 1.0, B, DIFFERENT_NONZERO_PATTERN));
 92:   PetscCall(MatSetOption(Bt, MAT_SYMMETRIC, PETSC_TRUE));
 93:   PetscCall(MatConvert(Bt, MATSBAIJ, MAT_INITIAL_MATRIX, &E));
 94:   PetscCall(MatDestroy(&E));
 95:   PetscCall(MatDestroy(&Bt));

 97:   /* Test MatIncreaseOverlap() */
 98:   PetscCall(PetscMalloc1(nd, &is1));
 99:   PetscCall(PetscMalloc1(nd, &is2));

101:   emptynd = PETSC_FALSE;
102:   if (rank == 0 && test_nd0) emptynd = PETSC_TRUE; /* test case */

104:   for (i = 0; i < nd; i++) {
105:     PetscCall(PetscRandomGetValue(rdm, &rval));
106:     sz = (int)(PetscRealPart(rval) * m);
107:     for (j = 0; j < sz; j++) {
108:       PetscCall(PetscRandomGetValue(rdm, &rval));
109:       idx[j * bs] = bs * (int)(PetscRealPart(rval) * Mbs);
110:       for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
111:     }
112:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is1 + i));
113:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is2 + i));
114:   }
115:   PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
116:   PetscCall(MatIncreaseOverlap(B, nd, is2, ov));

118:   for (i = 0; i < nd; ++i) {
119:     PetscCall(ISEqual(is1[i], is2[i], &flg));

121:     if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "i=%" PetscInt_FMT ", flg=%d :bs=%" PetscInt_FMT " m=%" PetscInt_FMT " ov=%" PetscInt_FMT " nd=%" PetscInt_FMT " np=%d\n", i, flg, bs, m, ov, nd, size));
122:   }

124:   for (i = 0; i < nd; ++i) {
125:     PetscCall(ISSort(is1[i]));
126:     PetscCall(ISSort(is2[i]));
127:   }

129:   PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB));
130:   PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));

132:   /* Test MatMult() */
133:   for (i = 0; i < nd; i++) {
134:     PetscCall(MatGetSize(submatA[i], &mm, &nn));
135:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
136:     PetscCall(VecDuplicate(xx, &s1));
137:     PetscCall(VecDuplicate(xx, &s2));
138:     for (j = 0; j < 3; j++) {
139:       PetscCall(VecSetRandom(xx, rdm));
140:       PetscCall(MatMult(submatA[i], xx, s1));
141:       PetscCall(MatMult(submatB[i], xx, s2));
142:       PetscCall(VecNorm(s1, NORM_2, &s1norm));
143:       PetscCall(VecNorm(s2, NORM_2, &s2norm));
144:       rnorm = s2norm - s1norm;
145:       if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm));
146:     }
147:     PetscCall(VecDestroy(&xx));
148:     PetscCall(VecDestroy(&s1));
149:     PetscCall(VecDestroy(&s2));
150:   }

152:   /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
153:   PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
154:   PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB));

156:   /* Test MatMult() */
157:   for (i = 0; i < nd; i++) {
158:     PetscCall(MatGetSize(submatA[i], &mm, &nn));
159:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
160:     PetscCall(VecDuplicate(xx, &s1));
161:     PetscCall(VecDuplicate(xx, &s2));
162:     for (j = 0; j < 3; j++) {
163:       PetscCall(VecSetRandom(xx, rdm));
164:       PetscCall(MatMult(submatA[i], xx, s1));
165:       PetscCall(MatMult(submatB[i], xx, s2));
166:       PetscCall(VecNorm(s1, NORM_2, &s1norm));
167:       PetscCall(VecNorm(s2, NORM_2, &s2norm));
168:       rnorm = s2norm - s1norm;
169:       if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm));
170:     }
171:     PetscCall(VecDestroy(&xx));
172:     PetscCall(VecDestroy(&s1));
173:     PetscCall(VecDestroy(&s2));
174:   }

176:   /* Free allocated memory */
177:   for (i = 0; i < nd; ++i) {
178:     PetscCall(ISDestroy(&is1[i]));
179:     PetscCall(ISDestroy(&is2[i]));
180:   }
181:   PetscCall(MatDestroySubMatrices(nd, &submatA));
182:   PetscCall(MatDestroySubMatrices(nd, &submatB));

184:   PetscCall(PetscFree(is1));
185:   PetscCall(PetscFree(is2));
186:   PetscCall(PetscFree(idx));
187:   PetscCall(PetscFree(rows));
188:   PetscCall(PetscFree(cols));
189:   PetscCall(PetscFree(vals));
190:   PetscCall(MatDestroy(&A));
191:   PetscCall(MatDestroy(&B));
192:   PetscCall(PetscRandomDestroy(&rdm));
193:   PetscCall(PetscFinalize());
194:   return 0;
195: }

197: /*TEST

199:    test:
200:       nsize: {{1 3}}
201:       args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd 7
202:       output_file: output/ex54.out

204:    test:
205:       suffix: 2
206:       args: -nd 2 -test_nd0
207:       output_file: output/ex54.out

209:    test:
210:       suffix: 3
211:       nsize: 3
212:       args: -nd 2 -test_nd0
213:       output_file: output/ex54.out

215: TEST*/