Actual source code: ex75.c

  1: static char help[] = "Tests various routines in MatMPISBAIJ format.\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec         x, y, u, s1, s2;
  8:   Mat         A, sA, sB;
  9:   PetscRandom rctx;
 10:   PetscReal   r1, r2, rnorm, tol = PETSC_SQRT_MACHINE_EPSILON;
 11:   PetscScalar one = 1.0, neg_one = -1.0, value[3], four = 4.0, alpha = 0.1;
 12:   PetscInt    n, col[3], n1, block, row, i, j, i2, j2, Ii, J, rstart, rend, bs = 1, mbs = 16, d_nz = 3, o_nz = 3, prob = 1;
 13:   PetscMPIInt size, rank;
 14:   PetscBool   flg;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 18:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
 19:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
 20:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-prob", &prob, NULL));

 22:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 23:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 25:   /* Create a BAIJ matrix A */
 26:   n = mbs * bs;
 27:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 28:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 29:   PetscCall(MatSetType(A, MATBAIJ));
 30:   PetscCall(MatSetFromOptions(A));
 31:   PetscCall(MatMPIBAIJSetPreallocation(A, bs, d_nz, NULL, o_nz, NULL));
 32:   PetscCall(MatSeqBAIJSetPreallocation(A, bs, d_nz, NULL));
 33:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));

 35:   if (bs == 1) {
 36:     if (prob == 1) { /* tridiagonal matrix */
 37:       value[0] = -1.0;
 38:       value[1] = 0.0;
 39:       value[2] = -1.0;
 40:       for (i = 1; i < n - 1; i++) {
 41:         col[0] = i - 1;
 42:         col[1] = i;
 43:         col[2] = i + 1;
 44:         PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 45:       }
 46:       i        = n - 1;
 47:       col[0]   = 0;
 48:       col[1]   = n - 2;
 49:       col[2]   = n - 1;
 50:       value[0] = 0.1;
 51:       value[1] = -1.0;
 52:       value[2] = 0.0;
 53:       PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));

 55:       i        = 0;
 56:       col[0]   = 0;
 57:       col[1]   = 1;
 58:       col[2]   = n - 1;
 59:       value[0] = 0.0;
 60:       value[1] = -1.0;
 61:       value[2] = 0.1;
 62:       PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 63:     } else if (prob == 2) { /* matrix for the five point stencil */
 64:       n1 = (int)PetscSqrtReal((PetscReal)n);
 65:       for (i = 0; i < n1; i++) {
 66:         for (j = 0; j < n1; j++) {
 67:           Ii = j + n1 * i;
 68:           if (i > 0) {
 69:             J = Ii - n1;
 70:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 71:           }
 72:           if (i < n1 - 1) {
 73:             J = Ii + n1;
 74:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 75:           }
 76:           if (j > 0) {
 77:             J = Ii - 1;
 78:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 79:           }
 80:           if (j < n1 - 1) {
 81:             J = Ii + 1;
 82:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 83:           }
 84:           PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES));
 85:         }
 86:       }
 87:     }
 88:     /* end of if (bs == 1) */
 89:   } else { /* bs > 1 */
 90:     for (block = 0; block < n / bs; block++) {
 91:       /* diagonal blocks */
 92:       value[0] = -1.0;
 93:       value[1] = 4.0;
 94:       value[2] = -1.0;
 95:       for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
 96:         col[0] = i - 1;
 97:         col[1] = i;
 98:         col[2] = i + 1;
 99:         PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
100:       }
101:       i        = bs - 1 + block * bs;
102:       col[0]   = bs - 2 + block * bs;
103:       col[1]   = bs - 1 + block * bs;
104:       value[0] = -1.0;
105:       value[1] = 4.0;
106:       PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));

108:       i        = 0 + block * bs;
109:       col[0]   = 0 + block * bs;
110:       col[1]   = 1 + block * bs;
111:       value[0] = 4.0;
112:       value[1] = -1.0;
113:       PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
114:     }
115:     /* off-diagonal blocks */
116:     value[0] = -1.0;
117:     for (i = 0; i < (n / bs - 1) * bs; i++) {
118:       col[0] = i + bs;
119:       PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
120:       col[0] = i;
121:       row    = i + bs;
122:       PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
123:     }
124:   }
125:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
126:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
127:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));

129:   /* Get SBAIJ matrix sA from A */
130:   PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));

132:   /* Test MatGetSize(), MatGetLocalSize() */
133:   PetscCall(MatGetSize(sA, &i, &j));
134:   PetscCall(MatGetSize(A, &i2, &j2));
135:   i -= i2;
136:   j -= j2;
137:   if (i || j) {
138:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetSize()\n", rank));
139:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
140:   }

142:   PetscCall(MatGetLocalSize(sA, &i, &j));
143:   PetscCall(MatGetLocalSize(A, &i2, &j2));
144:   i2 -= i;
145:   j2 -= j;
146:   if (i2 || j2) {
147:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetLocalSize()\n", rank));
148:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
149:   }

151:   /* vectors */
152:   /*--------------------*/
153:   /* i is obtained from MatGetLocalSize() */
154:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
155:   PetscCall(VecSetSizes(x, i, PETSC_DECIDE));
156:   PetscCall(VecSetFromOptions(x));
157:   PetscCall(VecDuplicate(x, &y));
158:   PetscCall(VecDuplicate(x, &u));
159:   PetscCall(VecDuplicate(x, &s1));
160:   PetscCall(VecDuplicate(x, &s2));

162:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
163:   PetscCall(PetscRandomSetFromOptions(rctx));
164:   PetscCall(VecSetRandom(x, rctx));
165:   PetscCall(VecSet(u, one));

167:   /* Test MatNorm() */
168:   PetscCall(MatNorm(A, NORM_FROBENIUS, &r1));
169:   PetscCall(MatNorm(sA, NORM_FROBENIUS, &r2));
170:   rnorm = PetscAbsReal(r1 - r2) / r2;
171:   if (rnorm > tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs));
172:   PetscCall(MatNorm(A, NORM_INFINITY, &r1));
173:   PetscCall(MatNorm(sA, NORM_INFINITY, &r2));
174:   rnorm = PetscAbsReal(r1 - r2) / r2;
175:   if (rnorm > tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs));
176:   PetscCall(MatNorm(A, NORM_1, &r1));
177:   PetscCall(MatNorm(sA, NORM_1, &r2));
178:   rnorm = PetscAbsReal(r1 - r2) / r2;
179:   if (rnorm > tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs));

181:   /* Test MatGetOwnershipRange() */
182:   PetscCall(MatGetOwnershipRange(sA, &rstart, &rend));
183:   PetscCall(MatGetOwnershipRange(A, &i2, &j2));
184:   i2 -= rstart;
185:   j2 -= rend;
186:   if (i2 || j2) {
187:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MaGetOwnershipRange()\n", rank));
188:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
189:   }

191:   /* Test MatDiagonalScale() */
192:   PetscCall(MatDiagonalScale(A, x, x));
193:   PetscCall(MatDiagonalScale(sA, x, x));
194:   PetscCall(MatMultEqual(A, sA, 10, &flg));
195:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDiagonalScale");

197:   /* Test MatGetDiagonal(), MatScale() */
198:   PetscCall(MatGetDiagonal(A, s1));
199:   PetscCall(MatGetDiagonal(sA, s2));
200:   PetscCall(VecNorm(s1, NORM_1, &r1));
201:   PetscCall(VecNorm(s2, NORM_1, &r2));
202:   r1 -= r2;
203:   if (r1 < -tol || r1 > tol) {
204:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n", rank, (double)r1));
205:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
206:   }

208:   PetscCall(MatScale(A, alpha));
209:   PetscCall(MatScale(sA, alpha));

211:   /* Test MatGetRowMaxAbs() */
212:   PetscCall(MatGetRowMaxAbs(A, s1, NULL));
213:   PetscCall(MatGetRowMaxAbs(sA, s2, NULL));

215:   PetscCall(VecNorm(s1, NORM_1, &r1));
216:   PetscCall(VecNorm(s2, NORM_1, &r2));
217:   r1 -= r2;
218:   if (r1 < -tol || r1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetRowMaxAbs() \n"));

220:   /* Test MatMult(), MatMultAdd() */
221:   PetscCall(MatMultEqual(A, sA, 10, &flg));
222:   if (!flg) {
223:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale()\n", rank));
224:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
225:   }

227:   PetscCall(MatMultAddEqual(A, sA, 10, &flg));
228:   if (!flg) {
229:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd()\n", rank));
230:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
231:   }

233:   /* Test MatMultTranspose(), MatMultTransposeAdd() */
234:   for (i = 0; i < 10; i++) {
235:     PetscCall(VecSetRandom(x, rctx));
236:     PetscCall(MatMultTranspose(A, x, s1));
237:     PetscCall(MatMultTranspose(sA, x, s2));
238:     PetscCall(VecNorm(s1, NORM_1, &r1));
239:     PetscCall(VecNorm(s2, NORM_1, &r2));
240:     r1 -= r2;
241:     if (r1 < -tol || r1 > tol) {
242:       PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale(), err=%g\n", rank, (double)r1));
243:       PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
244:     }
245:   }
246:   for (i = 0; i < 10; i++) {
247:     PetscCall(VecSetRandom(x, rctx));
248:     PetscCall(VecSetRandom(y, rctx));
249:     PetscCall(MatMultTransposeAdd(A, x, y, s1));
250:     PetscCall(MatMultTransposeAdd(sA, x, y, s2));
251:     PetscCall(VecNorm(s1, NORM_1, &r1));
252:     PetscCall(VecNorm(s2, NORM_1, &r2));
253:     r1 -= r2;
254:     if (r1 < -tol || r1 > tol) {
255:       PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd(), err=%g \n", rank, (double)r1));
256:       PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
257:     }
258:   }

260:   /* Test MatDuplicate() */
261:   PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &sB));
262:   PetscCall(MatEqual(sA, sB, &flg));
263:   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error in MatDuplicate(), sA != sB \n"));
264:   PetscCall(MatMultEqual(sA, sB, 5, &flg));
265:   if (!flg) {
266:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMult()\n", rank));
267:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
268:   }
269:   PetscCall(MatMultAddEqual(sA, sB, 5, &flg));
270:   if (!flg) {
271:     PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMultAdd(()\n", rank));
272:     PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
273:   }
274:   PetscCall(MatDestroy(&sB));
275:   PetscCall(VecDestroy(&u));
276:   PetscCall(VecDestroy(&x));
277:   PetscCall(VecDestroy(&y));
278:   PetscCall(VecDestroy(&s1));
279:   PetscCall(VecDestroy(&s2));
280:   PetscCall(MatDestroy(&sA));
281:   PetscCall(MatDestroy(&A));
282:   PetscCall(PetscRandomDestroy(&rctx));
283:   PetscCall(PetscFinalize());
284:   return 0;
285: }

287: /*TEST

289:    test:
290:       nsize: {{1 3}}
291:       args: -bs {{1 2 3  5  7 8}} -mat_ignore_lower_triangular -prob {{1 2}}

293: TEST*/