Actual source code: ex91.c


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

  4: #include <petscmat.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

205: /*TEST

207:    test:
208:       args: -ov 2

210: TEST*/