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;

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

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

 26:   /* Create a BAIJ matrix A */
 27:   n = mbs * bs;
 28:   MatCreate(PETSC_COMM_WORLD, &A);
 29:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
 30:   MatSetType(A, MATBAIJ);
 31:   MatSetFromOptions(A);
 32:   MatMPIBAIJSetPreallocation(A, bs, d_nz, NULL, o_nz, NULL);
 33:   MatSeqBAIJSetPreallocation(A, bs, d_nz, NULL);
 34:   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:         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:       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:       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:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 72:           }
 73:           if (i < n1 - 1) {
 74:             J = Ii + n1;
 75:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 76:           }
 77:           if (j > 0) {
 78:             J = Ii - 1;
 79:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 80:           }
 81:           if (j < n1 - 1) {
 82:             J = Ii + 1;
 83:             MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
 84:           }
 85:           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:         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:       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:       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:       MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES);
121:       col[0] = i;
122:       row    = i + bs;
123:       MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES);
124:     }
125:   }
126:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
127:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
128:   MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);

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

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

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

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

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

168:   /* Test MatNorm() */
169:   MatNorm(A, NORM_FROBENIUS, &r1);
170:   MatNorm(sA, NORM_FROBENIUS, &r2);
171:   rnorm = PetscAbsReal(r1 - r2) / r2;
172:   if (rnorm > tol && rank == 0) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs);
173:   MatNorm(A, NORM_INFINITY, &r1);
174:   MatNorm(sA, NORM_INFINITY, &r2);
175:   rnorm = PetscAbsReal(r1 - r2) / r2;
176:   if (rnorm > tol && rank == 0) PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs);
177:   MatNorm(A, NORM_1, &r1);
178:   MatNorm(sA, NORM_1, &r2);
179:   rnorm = PetscAbsReal(r1 - r2) / r2;
180:   if (rnorm > tol && rank == 0) 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:   MatGetOwnershipRange(sA, &rstart, &rend);
184:   MatGetOwnershipRange(A, &i2, &j2);
185:   i2 -= rstart;
186:   j2 -= rend;
187:   if (i2 || j2) {
188:     PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MaGetOwnershipRange()\n", rank);
189:     PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
190:   }

192:   /* Test MatDiagonalScale() */
193:   MatDiagonalScale(A, x, x);
194:   MatDiagonalScale(sA, x, x);
195:   MatMultEqual(A, sA, 10, &flg);

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

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

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

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

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

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

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

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