Actual source code: ex91.c

  1: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat          A, Atrans, sA, *submatA, *submatsA;
  8:   PetscInt     bs = 1, m = 43, ov = 1, i, j, k, *rows, *cols, M, nd = 5, *idx, mm, nn;
  9:   PetscMPIInt  size;
 10:   PetscScalar *vals, rval, one = 1.0;
 11:   IS          *is1, *is2;
 12:   PetscRandom  rand;
 13:   Vec          xx, s1, s2;
 14:   PetscReal    s1norm, s2norm, rnorm, tol = 10 * PETSC_SMALL;
 15:   PetscBool    flg;

 17:   PetscFunctionBeginUser;
 18:   PetscCall(PetscInitialize(&argc, &args, NULL, 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));

 24:   /* create a SeqBAIJ matrix A */
 25:   M = m * bs;
 26:   PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A));
 27:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
 28:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand));
 29:   PetscCall(PetscRandomSetFromOptions(rand));

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

 36:   /* Now set blocks of random values */
 37:   /* first, set diagonal blocks as zero */
 38:   for (j = 0; j < bs * bs; j++) vals[j] = 0.0;
 39:   for (i = 0; i < m; i++) {
 40:     cols[0] = i * bs;
 41:     rows[0] = i * bs;
 42:     for (j = 1; j < bs; j++) {
 43:       rows[j] = rows[j - 1] + 1;
 44:       cols[j] = cols[j - 1] + 1;
 45:     }
 46:     PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES));
 47:   }
 48:   /* second, add random blocks */
 49:   for (i = 0; i < 20 * bs; i++) {
 50:     PetscCall(PetscRandomGetValue(rand, &rval));
 51:     cols[0] = bs * (int)(PetscRealPart(rval) * m);
 52:     PetscCall(PetscRandomGetValue(rand, &rval));
 53:     rows[0] = bs * (int)(PetscRealPart(rval) * m);
 54:     for (j = 1; j < bs; j++) {
 55:       rows[j] = rows[j - 1] + 1;
 56:       cols[j] = cols[j - 1] + 1;
 57:     }

 59:     for (j = 0; j < bs * bs; j++) {
 60:       PetscCall(PetscRandomGetValue(rand, &rval));
 61:       vals[j] = rval;
 62:     }
 63:     PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES));
 64:   }

 66:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 67:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 69:   /* make A a symmetric matrix: A <- A^T + A */
 70:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
 71:   PetscCall(MatAXPY(A, one, Atrans, DIFFERENT_NONZERO_PATTERN));
 72:   PetscCall(MatDestroy(&Atrans));
 73:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
 74:   PetscCall(MatEqual(A, Atrans, &flg));
 75:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A+A^T is non-symmetric");
 76:   PetscCall(MatDestroy(&Atrans));

 78:   /* create a SeqSBAIJ matrix sA (= A) */
 79:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
 80:   PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA));

 82:   /* Test sA==A through MatMult() */
 83:   for (i = 0; i < nd; i++) {
 84:     PetscCall(MatGetSize(A, &mm, &nn));
 85:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
 86:     PetscCall(VecDuplicate(xx, &s1));
 87:     PetscCall(VecDuplicate(xx, &s2));
 88:     for (j = 0; j < 3; j++) {
 89:       PetscCall(VecSetRandom(xx, rand));
 90:       PetscCall(MatMult(A, xx, s1));
 91:       PetscCall(MatMult(sA, xx, s2));
 92:       PetscCall(VecNorm(s1, NORM_2, &s1norm));
 93:       PetscCall(VecNorm(s2, NORM_2, &s2norm));
 94:       rnorm = s2norm - s1norm;
 95:       if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm));
 96:     }
 97:     PetscCall(VecDestroy(&xx));
 98:     PetscCall(VecDestroy(&s1));
 99:     PetscCall(VecDestroy(&s2));
100:   }

102:   /* Test MatIncreaseOverlap() */
103:   PetscCall(PetscMalloc1(nd, &is1));
104:   PetscCall(PetscMalloc1(nd, &is2));

106:   for (i = 0; i < nd; i++) {
107:     PetscCall(PetscRandomGetValue(rand, &rval));
108:     size = (int)(PetscRealPart(rval) * m);
109:     for (j = 0; j < size; j++) {
110:       PetscCall(PetscRandomGetValue(rand, &rval));
111:       idx[j * bs] = bs * (int)(PetscRealPart(rval) * m);
112:       for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
113:     }
114:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size * bs, idx, PETSC_COPY_VALUES, is1 + i));
115:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size * bs, idx, PETSC_COPY_VALUES, is2 + i));
116:   }
117:   /* for debugging */
118:   /*
119:   PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF));
120:   PetscCall(MatView(sA,PETSC_VIEWER_STDOUT_SELF));
121:   */

123:   PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
124:   PetscCall(MatIncreaseOverlap(sA, nd, is2, ov));

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

131:   for (i = 0; i < nd; ++i) {
132:     PetscCall(ISEqual(is1[i], is2[i], &flg));
133:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", is1 != is2", i);
134:   }

136:   PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
137:   PetscCall(MatCreateSubMatrices(sA, nd, is2, is2, MAT_INITIAL_MATRIX, &submatsA));

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

159:   /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
160:   PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
161:   PetscCall(MatCreateSubMatrices(sA, nd, is2, is2, MAT_REUSE_MATRIX, &submatsA));

163:   /* Test MatMult() */
164:   for (i = 0; i < nd; i++) {
165:     PetscCall(MatGetSize(submatA[i], &mm, &nn));
166:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
167:     PetscCall(VecDuplicate(xx, &s1));
168:     PetscCall(VecDuplicate(xx, &s2));
169:     for (j = 0; j < 3; j++) {
170:       PetscCall(VecSetRandom(xx, rand));
171:       PetscCall(MatMult(submatA[i], xx, s1));
172:       PetscCall(MatMult(submatsA[i], xx, s2));
173:       PetscCall(VecNorm(s1, NORM_2, &s1norm));
174:       PetscCall(VecNorm(s2, NORM_2, &s2norm));
175:       rnorm = s2norm - s1norm;
176:       if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm));
177:     }
178:     PetscCall(VecDestroy(&xx));
179:     PetscCall(VecDestroy(&s1));
180:     PetscCall(VecDestroy(&s2));
181:   }

183:   /* Free allocated memory */
184:   for (i = 0; i < nd; ++i) {
185:     PetscCall(ISDestroy(&is1[i]));
186:     PetscCall(ISDestroy(&is2[i]));
187:   }
188:   PetscCall(MatDestroySubMatrices(nd, &submatA));
189:   PetscCall(MatDestroySubMatrices(nd, &submatsA));

191:   PetscCall(PetscFree(is1));
192:   PetscCall(PetscFree(is2));
193:   PetscCall(PetscFree(idx));
194:   PetscCall(PetscFree(rows));
195:   PetscCall(PetscFree(cols));
196:   PetscCall(PetscFree(vals));
197:   PetscCall(MatDestroy(&A));
198:   PetscCall(MatDestroy(&sA));
199:   PetscCall(PetscRandomDestroy(&rand));
200:   PetscCall(PetscFinalize());
201:   return 0;
202: }

204: /*TEST

206:    test:
207:       args: -ov 2

209: TEST*/