Actual source code: ex77.c


  2: static char help[] = "Tests the various sequential routines in MatSBAIJ format. Same as ex74.c except diagonal entries of the matrices are zeros.\n";

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Vec             x, y, b, s1, s2;
  9:   Mat             A;  /* linear system matrix */
 10:   Mat             sA; /* symmetric part of the matrices */
 11:   PetscInt        n, mbs = 16, bs = 1, nz = 3, prob = 2, i, j, col[3], row, Ii, J, n1;
 12:   const PetscInt *ip_ptr;
 13:   PetscScalar     neg_one = -1.0, value[3], alpha = 0.1;
 14:   PetscMPIInt     size;
 15:   IS              ip, isrow, iscol;
 16:   PetscRandom     rdm;
 17:   PetscBool       reorder = PETSC_FALSE;
 18:   MatInfo         minfo1, minfo2;
 19:   PetscReal       norm1, norm2, tol = 10 * PETSC_SMALL;

 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 23:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 24:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 25:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
 26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));

 28:   n = mbs * bs;
 29:   PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &A));
 30:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
 31:   PetscCall(MatCreateSeqSBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &sA));
 32:   PetscCall(MatSetOption(sA, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));

 34:   /* Test MatGetOwnershipRange() */
 35:   PetscCall(MatGetOwnershipRange(A, &Ii, &J));
 36:   PetscCall(MatGetOwnershipRange(sA, &i, &j));
 37:   if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetOwnershipRange() in MatSBAIJ format\n"));

 39:   /* Assemble matrix */
 40:   if (bs == 1) {
 41:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL));
 42:     if (prob == 1) { /* tridiagonal matrix */
 43:       value[0] = -1.0;
 44:       value[1] = 2.0;
 45:       value[2] = -1.0;
 46:       for (i = 1; i < n - 1; i++) {
 47:         col[0] = i - 1;
 48:         col[1] = i;
 49:         col[2] = i + 1;
 50:         PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 51:         PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
 52:       }
 53:       i      = n - 1;
 54:       col[0] = 0;
 55:       col[1] = n - 2;
 56:       col[2] = n - 1;

 58:       value[0] = 0.1;
 59:       value[1] = -1;
 60:       value[2] = 2;
 61:       PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 62:       PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));

 64:       i      = 0;
 65:       col[0] = 0;
 66:       col[1] = 1;
 67:       col[2] = n - 1;

 69:       value[0] = 2.0;
 70:       value[1] = -1.0;
 71:       value[2] = 0.1;
 72:       PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 73:       PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
 74:     } else if (prob == 2) { /* matrix for the five point stencil */
 75:       n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001);
 76:       PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!");
 77:       for (i = 0; i < n1; i++) {
 78:         for (j = 0; j < n1; j++) {
 79:           Ii = j + n1 * i;
 80:           if (i > 0) {
 81:             J = Ii - n1;
 82:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 83:             PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 84:           }
 85:           if (i < n1 - 1) {
 86:             J = Ii + n1;
 87:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 88:             PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 89:           }
 90:           if (j > 0) {
 91:             J = Ii - 1;
 92:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 93:             PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 94:           }
 95:           if (j < n1 - 1) {
 96:             J = Ii + 1;
 97:             PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 98:             PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
 99:           }
100:         }
101:       }
102:     }
103:   } else { /* bs > 1 */
104: #if defined(DIAGB)
105:     for (block = 0; block < n / bs; block++) {
106:       /* diagonal blocks */
107:       value[0] = -1.0;
108:       value[1] = 4.0;
109:       value[2] = -1.0;
110:       for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
111:         col[0] = i - 1;
112:         col[1] = i;
113:         col[2] = i + 1;
114:         PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
115:         PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
116:       }
117:       i      = bs - 1 + block * bs;
118:       col[0] = bs - 2 + block * bs;
119:       col[1] = bs - 1 + block * bs;

121:       value[0] = -1.0;
122:       value[1] = 4.0;
123:       PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
124:       PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES));

126:       i      = 0 + block * bs;
127:       col[0] = 0 + block * bs;
128:       col[1] = 1 + block * bs;

130:       value[0] = 4.0;
131:       value[1] = -1.0;
132:       PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
133:       PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES));
134:     }
135: #endif
136:     /* off-diagonal blocks */
137:     value[0] = -1.0;
138:     for (i = 0; i < (n / bs - 1) * bs; i++) {
139:       col[0] = i + bs;
140:       PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
141:       PetscCall(MatSetValues(sA, 1, &i, 1, col, value, INSERT_VALUES));
142:       col[0] = i;
143:       row    = i + bs;
144:       PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
145:       PetscCall(MatSetValues(sA, 1, &row, 1, col, value, INSERT_VALUES));
146:     }
147:   }
148:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
149:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

151:   PetscCall(MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY));
152:   PetscCall(MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY));

154:   /* Test MatNorm() */
155:   PetscCall(MatNorm(A, NORM_FROBENIUS, &norm1));
156:   PetscCall(MatNorm(sA, NORM_FROBENIUS, &norm2));
157:   norm1 -= norm2;
158:   if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm(), fnorm1-fnorm2=%16.14e\n", (double)norm1));
159:   PetscCall(MatNorm(A, NORM_INFINITY, &norm1));
160:   PetscCall(MatNorm(sA, NORM_INFINITY, &norm2));
161:   norm1 -= norm2;
162:   if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n", (double)norm1));

164:   /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
165:   PetscCall(MatGetInfo(A, MAT_LOCAL, &minfo1));
166:   PetscCall(MatGetInfo(sA, MAT_LOCAL, &minfo2));
167:   i = (int)(minfo1.nz_used - minfo2.nz_used);
168:   j = (int)(minfo1.nz_allocated - minfo2.nz_allocated);
169:   if (i < 0 || j < 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetInfo()\n"));

171:   PetscCall(MatGetSize(A, &Ii, &J));
172:   PetscCall(MatGetSize(sA, &i, &j));
173:   if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetSize()\n"));

175:   PetscCall(MatGetBlockSize(A, &Ii));
176:   PetscCall(MatGetBlockSize(sA, &i));
177:   if (i - Ii) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetBlockSize()\n"));

179:   /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
180:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
181:   PetscCall(PetscRandomSetFromOptions(rdm));
182:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
183:   PetscCall(VecDuplicate(x, &s1));
184:   PetscCall(VecDuplicate(x, &s2));
185:   PetscCall(VecDuplicate(x, &y));
186:   PetscCall(VecDuplicate(x, &b));

188:   PetscCall(VecSetRandom(x, rdm));

190:   PetscCall(MatDiagonalScale(A, x, x));
191:   PetscCall(MatDiagonalScale(sA, x, x));

193:   PetscCall(MatGetDiagonal(A, s1));
194:   PetscCall(MatGetDiagonal(sA, s2));
195:   PetscCall(VecNorm(s1, NORM_1, &norm1));
196:   PetscCall(VecNorm(s2, NORM_1, &norm2));
197:   norm1 -= norm2;
198:   if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatGetDiagonal() \n"));

200:   PetscCall(MatScale(A, alpha));
201:   PetscCall(MatScale(sA, alpha));

203:   /* Test MatMult(), MatMultAdd() */
204:   for (i = 0; i < 40; i++) {
205:     PetscCall(VecSetRandom(x, rdm));
206:     PetscCall(MatMult(A, x, s1));
207:     PetscCall(MatMult(sA, x, s2));
208:     PetscCall(VecNorm(s1, NORM_1, &norm1));
209:     PetscCall(VecNorm(s2, NORM_1, &norm2));
210:     norm1 -= norm2;
211:     if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMult(), MatDiagonalScale() or MatScale()\n"));
212:   }

214:   for (i = 0; i < 40; i++) {
215:     PetscCall(VecSetRandom(x, rdm));
216:     PetscCall(VecSetRandom(y, rdm));
217:     PetscCall(MatMultAdd(A, x, y, s1));
218:     PetscCall(MatMultAdd(sA, x, y, s2));
219:     PetscCall(VecNorm(s1, NORM_1, &norm1));
220:     PetscCall(VecNorm(s2, NORM_1, &norm2));
221:     norm1 -= norm2;
222:     if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n"));
223:   }

225:   /* Test MatReordering() */
226:   PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &isrow, &iscol));
227:   ip = isrow;

229:   if (reorder) {
230:     IS        nip;
231:     PetscInt *nip_ptr;
232:     PetscCall(PetscMalloc1(mbs, &nip_ptr));
233:     PetscCall(ISGetIndices(ip, &ip_ptr));
234:     PetscCall(PetscArraycpy(nip_ptr, ip_ptr, mbs));
235:     i                = nip_ptr[1];
236:     nip_ptr[1]       = nip_ptr[mbs - 2];
237:     nip_ptr[mbs - 2] = i;
238:     i                = nip_ptr[0];
239:     nip_ptr[0]       = nip_ptr[mbs - 1];
240:     nip_ptr[mbs - 1] = i;
241:     PetscCall(ISRestoreIndices(ip, &ip_ptr));
242:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, mbs, nip_ptr, PETSC_COPY_VALUES, &nip));
243:     PetscCall(PetscFree(nip_ptr));

245:     PetscCall(MatReorderingSeqSBAIJ(sA, ip));
246:     PetscCall(ISDestroy(&nip));
247:   }

249:   PetscCall(ISDestroy(&iscol));
250:   PetscCall(ISDestroy(&isrow));
251:   PetscCall(MatDestroy(&A));
252:   PetscCall(MatDestroy(&sA));
253:   PetscCall(VecDestroy(&x));
254:   PetscCall(VecDestroy(&y));
255:   PetscCall(VecDestroy(&s1));
256:   PetscCall(VecDestroy(&s2));
257:   PetscCall(VecDestroy(&b));
258:   PetscCall(PetscRandomDestroy(&rdm));

260:   PetscCall(PetscFinalize());
261:   return 0;
262: }

264: /*TEST

266:    test:
267:       args: -bs {{1 2 3 4 5 6 7 8}}

269: TEST*/