Actual source code: ex75.c


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

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Vec         x, y, u, s1, s2;
  9:   Mat         A, sA, sB;
 10:   PetscRandom rctx;
 11:   PetscReal   r1, r2, rnorm, tol = PETSC_SQRT_MACHINE_EPSILON;
 12:   PetscScalar one = 1.0, neg_one = -1.0, value[3], four = 4.0, alpha = 0.1;
 13:   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;
 14:   PetscMPIInt size, rank;
 15:   PetscBool   flg;

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

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

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

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

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

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

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

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

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

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

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

168:   /* Test MatNorm() */
169:   PetscCall(MatNorm(A, NORM_FROBENIUS, &r1));
170:   PetscCall(MatNorm(sA, NORM_FROBENIUS, &r2));
171:   rnorm = PetscAbsReal(r1 - r2) / r2;
172:   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));
173:   PetscCall(MatNorm(A, NORM_INFINITY, &r1));
174:   PetscCall(MatNorm(sA, NORM_INFINITY, &r2));
175:   rnorm = PetscAbsReal(r1 - r2) / r2;
176:   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));
177:   PetscCall(MatNorm(A, NORM_1, &r1));
178:   PetscCall(MatNorm(sA, NORM_1, &r2));
179:   rnorm = PetscAbsReal(r1 - r2) / r2;
180:   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));

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

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

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

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

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

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

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

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

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

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

288: /*TEST

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

294: TEST*/