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;

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

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

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

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

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

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

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

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

 81:   for (i = 0; i < nd; i++) {
 82:     PetscRandomGetValue(rdm, &rval);
 83:     lsize = (int)(PetscRealPart(rval) * m);
 84:     for (j = 0; j < lsize; j++) {
 85:       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:     ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is1 + i);
 90:     ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is2 + i);
 91:   }
 92:   MatIncreaseOverlap(A, nd, is1, ov);
 93:   MatIncreaseOverlap(B, nd, is2, ov);

 95:   for (i = 0; i < nd; ++i) {
 96:     ISEqual(is1[i], is2[i], &flg);
 98:   }

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

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

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

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

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