Actual source code: ex54.c


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

  4: #include <petscmat.h>

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

 19:   PetscInitialize(&argc, &args, (char *)0, help);
 20:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 21:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

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

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

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

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

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

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

 59:   /* Now set blocks of values */
 60:   for (i = 0; i < 40 * bs; i++) {
 61:     PetscRandomGetValue(rdm, &rval);
 62:     cols[0] = bs * (int)(PetscRealPart(rval) * Mbs);
 63:     PetscRandomGetValue(rdm, &rval);
 64:     rows[0] = rstart + bs * (int)(PetscRealPart(rval) * m);
 65:     for (j = 1; j < bs; j++) {
 66:       rows[j] = rows[j - 1] + 1;
 67:       cols[j] = cols[j - 1] + 1;
 68:     }

 70:     for (j = 0; j < bs * bs; j++) {
 71:       PetscRandomGetValue(rdm, &rval);
 72:       vals[j] = rval;
 73:     }
 74:     MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES);
 75:     MatSetValues(B, bs, rows, bs, cols, vals, ADD_VALUES);
 76:   }
 77:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 78:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 79:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 80:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

 82:   /* Test MatIncreaseOverlap() */
 83:   PetscMalloc1(nd, &is1);
 84:   PetscMalloc1(nd, &is2);

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

 89:   for (i = 0; i < nd; i++) {
 90:     PetscRandomGetValue(rdm, &rval);
 91:     sz = (int)(PetscRealPart(rval) * m);
 92:     for (j = 0; j < sz; j++) {
 93:       PetscRandomGetValue(rdm, &rval);
 94:       idx[j * bs] = bs * (int)(PetscRealPart(rval) * Mbs);
 95:       for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
 96:     }
 97:     ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is1 + i);
 98:     ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is2 + i);
 99:   }
100:   MatIncreaseOverlap(A, nd, is1, ov);
101:   MatIncreaseOverlap(B, nd, is2, ov);

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

106:     if (!flg) 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);
107:   }

109:   for (i = 0; i < nd; ++i) {
110:     ISSort(is1[i]);
111:     ISSort(is2[i]);
112:   }

114:   MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB);
115:   MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA);

117:   /* Test MatMult() */
118:   for (i = 0; i < nd; i++) {
119:     MatGetSize(submatA[i], &mm, &nn);
120:     VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
121:     VecDuplicate(xx, &s1);
122:     VecDuplicate(xx, &s2);
123:     for (j = 0; j < 3; j++) {
124:       VecSetRandom(xx, rdm);
125:       MatMult(submatA[i], xx, s1);
126:       MatMult(submatB[i], xx, s2);
127:       VecNorm(s1, NORM_2, &s1norm);
128:       VecNorm(s2, NORM_2, &s2norm);
129:       rnorm = s2norm - s1norm;
130:       if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm);
131:     }
132:     VecDestroy(&xx);
133:     VecDestroy(&s1);
134:     VecDestroy(&s2);
135:   }

137:   /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
138:   MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA);
139:   MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB);

141:   /* Test MatMult() */
142:   for (i = 0; i < nd; i++) {
143:     MatGetSize(submatA[i], &mm, &nn);
144:     VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
145:     VecDuplicate(xx, &s1);
146:     VecDuplicate(xx, &s2);
147:     for (j = 0; j < 3; j++) {
148:       VecSetRandom(xx, rdm);
149:       MatMult(submatA[i], xx, s1);
150:       MatMult(submatB[i], xx, s2);
151:       VecNorm(s1, NORM_2, &s1norm);
152:       VecNorm(s2, NORM_2, &s2norm);
153:       rnorm = s2norm - s1norm;
154:       if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm);
155:     }
156:     VecDestroy(&xx);
157:     VecDestroy(&s1);
158:     VecDestroy(&s2);
159:   }

161:   /* Free allocated memory */
162:   for (i = 0; i < nd; ++i) {
163:     ISDestroy(&is1[i]);
164:     ISDestroy(&is2[i]);
165:   }
166:   MatDestroySubMatrices(nd, &submatA);
167:   MatDestroySubMatrices(nd, &submatB);

169:   PetscFree(is1);
170:   PetscFree(is2);
171:   PetscFree(idx);
172:   PetscFree(rows);
173:   PetscFree(cols);
174:   PetscFree(vals);
175:   MatDestroy(&A);
176:   MatDestroy(&B);
177:   PetscRandomDestroy(&rdm);
178:   PetscFinalize();
179:   return 0;
180: }

182: /*TEST

184:    test:
185:       nsize: {{1 3}}
186:       args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd 7
187:       output_file: output/ex54.out

189:    test:
190:       suffix: 2
191:       args: -nd 2 -test_nd0
192:       output_file: output/ex54.out

194:    test:
195:       suffix: 3
196:       nsize: 3
197:       args: -nd 2 -test_nd0
198:       output_file: output/ex54.out

200: TEST*/