Actual source code: ex77.c

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

  3: #include <petscmat.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

263: /*TEST

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

268: TEST*/