Actual source code: ex51.c


  2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for MatBAIJ format.\n";

  4: #include <petscmat.h>

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

 17:   PetscFunctionBeginUser;
 18:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 19:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
 20:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL));
 21:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
 22:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
 23:   M = m * bs;

 25:   PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A));
 26:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
 27:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, M, M, 15, NULL, &B));
 28:   PetscCall(MatSetBlockSize(B, bs));
 29:   PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
 30:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
 31:   PetscCall(PetscRandomSetFromOptions(rdm));

 33:   PetscCall(PetscMalloc1(bs, &rows));
 34:   PetscCall(PetscMalloc1(bs, &cols));
 35:   PetscCall(PetscMalloc1(bs * bs, &vals));
 36:   PetscCall(PetscMalloc1(M, &idx));

 38:   /* Now set blocks of values */
 39:   for (i = 0; i < 20 * bs; i++) {
 40:     PetscInt nr = 1, nc = 1;
 41:     PetscCall(PetscRandomGetValue(rdm, &rval));
 42:     cols[0] = bs * (int)(PetscRealPart(rval) * m);
 43:     PetscCall(PetscRandomGetValue(rdm, &rval));
 44:     rows[0] = bs * (int)(PetscRealPart(rval) * m);
 45:     for (j = 1; j < bs; j++) {
 46:       PetscCall(PetscRandomGetValue(rdm, &rval));
 47:       if (PetscRealPart(rval) > .5) rows[nr++] = rows[0] + j - 1;
 48:     }
 49:     for (j = 1; j < bs; j++) {
 50:       PetscCall(PetscRandomGetValue(rdm, &rval));
 51:       if (PetscRealPart(rval) > .5) cols[nc++] = cols[0] + j - 1;
 52:     }

 54:     for (j = 0; j < nr * nc; j++) {
 55:       PetscCall(PetscRandomGetValue(rdm, &rval));
 56:       vals[j] = rval;
 57:     }
 58:     PetscCall(MatSetValues(A, nr, rows, nc, cols, vals, ADD_VALUES));
 59:     PetscCall(MatSetValues(B, nr, rows, nc, cols, vals, ADD_VALUES));
 60:   }

 62:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 63:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 64:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 65:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

 67:   /* Test MatConvert_SeqAIJ_Seq(S)BAIJ handles incompletely filled blocks */
 68:   PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &E));
 69:   PetscCall(MatDestroy(&E));
 70:   PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &Bt));
 71:   PetscCall(MatAXPY(Bt, 1.0, B, DIFFERENT_NONZERO_PATTERN));
 72:   PetscCall(MatSetOption(Bt, MAT_SYMMETRIC, PETSC_TRUE));
 73:   PetscCall(MatConvert(Bt, MATSBAIJ, MAT_INITIAL_MATRIX, &E));
 74:   PetscCall(MatDestroy(&E));
 75:   PetscCall(MatDestroy(&Bt));

 77:   /* Test MatIncreaseOverlap() */
 78:   PetscCall(PetscMalloc1(nd, &is1));
 79:   PetscCall(PetscMalloc1(nd, &is2));

 81:   for (i = 0; i < nd; i++) {
 82:     PetscCall(PetscRandomGetValue(rdm, &rval));
 83:     lsize = (int)(PetscRealPart(rval) * m);
 84:     for (j = 0; j < lsize; j++) {
 85:       PetscCall(PetscRandomGetValue(rdm, &rval));
 86:       idx[j * bs] = bs * (int)(PetscRealPart(rval) * m);
 87:       for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
 88:     }
 89:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is1 + i));
 90:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is2 + i));
 91:   }
 92:   PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
 93:   PetscCall(MatIncreaseOverlap(B, nd, is2, ov));

 95:   for (i = 0; i < nd; ++i) {
 96:     PetscCall(ISEqual(is1[i], is2[i], &flg));
 97:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", flg =%d", i, (int)flg);
 98:   }

100:   for (i = 0; i < nd; ++i) {
101:     PetscCall(ISSort(is1[i]));
102:     PetscCall(ISSort(is2[i]));
103:   }

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

108:   /* Test MatMult() */
109:   for (i = 0; i < nd; i++) {
110:     PetscCall(MatGetSize(submatA[i], &mm, &nn));
111:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
112:     PetscCall(VecDuplicate(xx, &s1));
113:     PetscCall(VecDuplicate(xx, &s2));
114:     for (j = 0; j < 3; j++) {
115:       PetscCall(VecSetRandom(xx, rdm));
116:       PetscCall(MatMult(submatA[i], xx, s1));
117:       PetscCall(MatMult(submatB[i], xx, s2));
118:       PetscCall(VecNorm(s1, NORM_2, &s1norm));
119:       PetscCall(VecNorm(s2, NORM_2, &s2norm));
120:       rnorm = s2norm - s1norm;
121:       if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm));
122:     }
123:     PetscCall(VecDestroy(&xx));
124:     PetscCall(VecDestroy(&s1));
125:     PetscCall(VecDestroy(&s2));
126:   }
127:   /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
128:   PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
129:   PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB));

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

151:   /* Free allocated memory */
152:   for (i = 0; i < nd; ++i) {
153:     PetscCall(ISDestroy(&is1[i]));
154:     PetscCall(ISDestroy(&is2[i]));
155:   }
156:   PetscCall(MatDestroySubMatrices(nd, &submatA));
157:   PetscCall(MatDestroySubMatrices(nd, &submatB));
158:   PetscCall(PetscFree(is1));
159:   PetscCall(PetscFree(is2));
160:   PetscCall(PetscFree(idx));
161:   PetscCall(PetscFree(rows));
162:   PetscCall(PetscFree(cols));
163:   PetscCall(PetscFree(vals));
164:   PetscCall(MatDestroy(&A));
165:   PetscCall(MatDestroy(&B));
166:   PetscCall(PetscRandomDestroy(&rdm));
167:   PetscCall(PetscFinalize());
168:   return 0;
169: }

171: /*TEST

173:    test:
174:       args: -mat_block_size {{1 2  5 7 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}}

176: TEST*/