Actual source code: ex53.c

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

  3: #include <petscmat.h>
  4: #define IMAX 15
  5: int main(int argc, char **args)
  6: {
  7:   Mat                A, B, C, At, Bt;
  8:   PetscViewer        fd;
  9:   char               file[PETSC_MAX_PATH_LEN];
 10:   PetscRandom        rand;
 11:   Vec                xx, yy, s1, s2;
 12:   PetscReal          s1norm, s2norm, rnorm, tol = 1.e-10;
 13:   PetscInt           rstart, rend, rows[2], cols[2], m, n, i, j, M, N, ct, row, ncols1, ncols2, bs;
 14:   PetscMPIInt        rank, size;
 15:   const PetscInt    *cols1, *cols2;
 16:   PetscScalar        vals1[4], vals2[4], v;
 17:   const PetscScalar *v1, *v2;
 18:   PetscBool          flg;

 20:   PetscFunctionBeginUser;
 21:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 22:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 23:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 25:   /* Check out if MatLoad() works */
 26:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
 27:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Input file not specified");
 28:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 29:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 30:   PetscCall(MatSetType(A, MATBAIJ));
 31:   PetscCall(MatLoad(A, fd));

 33:   PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
 34:   PetscCall(PetscViewerDestroy(&fd));

 36:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 37:   PetscCall(PetscRandomSetFromOptions(rand));
 38:   PetscCall(MatGetLocalSize(A, &m, &n));
 39:   PetscCall(VecCreate(PETSC_COMM_WORLD, &xx));
 40:   PetscCall(VecSetSizes(xx, m, PETSC_DECIDE));
 41:   PetscCall(VecSetFromOptions(xx));
 42:   PetscCall(VecDuplicate(xx, &s1));
 43:   PetscCall(VecDuplicate(xx, &s2));
 44:   PetscCall(VecDuplicate(xx, &yy));
 45:   PetscCall(MatGetBlockSize(A, &bs));

 47:   /* Test MatNorm() */
 48:   PetscCall(MatNorm(A, NORM_FROBENIUS, &s1norm));
 49:   PetscCall(MatNorm(B, NORM_FROBENIUS, &s2norm));
 50:   rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
 51:   if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
 52:   PetscCall(MatNorm(A, NORM_INFINITY, &s1norm));
 53:   PetscCall(MatNorm(B, NORM_INFINITY, &s2norm));
 54:   rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
 55:   if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
 56:   PetscCall(MatNorm(A, NORM_1, &s1norm));
 57:   PetscCall(MatNorm(B, NORM_1, &s2norm));
 58:   rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm;
 59:   if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));

 61:   /* Test MatMult() */
 62:   for (i = 0; i < IMAX; i++) {
 63:     PetscCall(VecSetRandom(xx, rand));
 64:     PetscCall(MatMult(A, xx, s1));
 65:     PetscCall(MatMult(B, xx, s2));
 66:     PetscCall(VecAXPY(s2, -1.0, s1));
 67:     PetscCall(VecNorm(s2, NORM_2, &rnorm));
 68:     if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMult - Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)rnorm, bs));
 69:   }

 71:   /* test MatMultAdd() */
 72:   for (i = 0; i < IMAX; i++) {
 73:     PetscCall(VecSetRandom(xx, rand));
 74:     PetscCall(VecSetRandom(yy, rand));
 75:     PetscCall(MatMultAdd(A, xx, yy, s1));
 76:     PetscCall(MatMultAdd(B, xx, yy, s2));
 77:     PetscCall(VecAXPY(s2, -1.0, s1));
 78:     PetscCall(VecNorm(s2, NORM_2, &rnorm));
 79:     if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultAdd - Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)rnorm, bs));
 80:   }

 82:   /* Test MatMultTranspose() */
 83:   for (i = 0; i < IMAX; i++) {
 84:     PetscCall(VecSetRandom(xx, rand));
 85:     PetscCall(MatMultTranspose(A, xx, s1));
 86:     PetscCall(MatMultTranspose(B, xx, s2));
 87:     PetscCall(VecNorm(s1, NORM_2, &s1norm));
 88:     PetscCall(VecNorm(s2, NORM_2, &s2norm));
 89:     rnorm = s2norm - s1norm;
 90:     if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
 91:   }
 92:   /* Test MatMultTransposeAdd() */
 93:   for (i = 0; i < IMAX; i++) {
 94:     PetscCall(VecSetRandom(xx, rand));
 95:     PetscCall(VecSetRandom(yy, rand));
 96:     PetscCall(MatMultTransposeAdd(A, xx, yy, s1));
 97:     PetscCall(MatMultTransposeAdd(B, xx, yy, s2));
 98:     PetscCall(VecNorm(s1, NORM_2, &s1norm));
 99:     PetscCall(VecNorm(s2, NORM_2, &s2norm));
100:     rnorm = s2norm - s1norm;
101:     if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
102:   }

104:   /* Check MatGetValues() */
105:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
106:   PetscCall(MatGetSize(A, &M, &N));

108:   for (i = 0; i < IMAX; i++) {
109:     /* Create random row numbers ad col numbers */
110:     PetscCall(PetscRandomGetValue(rand, &v));
111:     cols[0] = (int)(PetscRealPart(v) * N);
112:     PetscCall(PetscRandomGetValue(rand, &v));
113:     cols[1] = (int)(PetscRealPart(v) * N);
114:     PetscCall(PetscRandomGetValue(rand, &v));
115:     rows[0] = rstart + (int)(PetscRealPart(v) * m);
116:     PetscCall(PetscRandomGetValue(rand, &v));
117:     rows[1] = rstart + (int)(PetscRealPart(v) * m);

119:     PetscCall(MatGetValues(A, 2, rows, 2, cols, vals1));
120:     PetscCall(MatGetValues(B, 2, rows, 2, cols, vals2));

122:     for (j = 0; j < 4; j++) {
123:       if (vals1[j] != vals2[j]) {
124:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]: Error: MatGetValues rstart = %2" PetscInt_FMT "  row = %2" PetscInt_FMT " col = %2" PetscInt_FMT " val1 = %e val2 = %e bs = %" PetscInt_FMT "\n", rank, rstart, rows[j / 2], cols[j % 2], (double)PetscRealPart(vals1[j]), (double)PetscRealPart(vals2[j]), bs));
125:       }
126:     }
127:   }

129:   /* Test MatGetRow()/ MatRestoreRow() */
130:   for (ct = 0; ct < 100; ct++) {
131:     PetscCall(PetscRandomGetValue(rand, &v));
132:     row = rstart + (PetscInt)(PetscRealPart(v) * m);
133:     PetscCall(MatGetRow(A, row, &ncols1, &cols1, &v1));
134:     PetscCall(MatGetRow(B, row, &ncols2, &cols2, &v2));

136:     for (i = 0, j = 0; i < ncols1 && j < ncols2; j++) {
137:       while (cols2[j] != cols1[i]) i++;
138:       PetscCheck(v1[i] == v2[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - vals incorrect.");
139:     }
140:     PetscCheck(j >= ncols2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - cols incorrect");

142:     PetscCall(MatRestoreRow(A, row, &ncols1, &cols1, &v1));
143:     PetscCall(MatRestoreRow(B, row, &ncols2, &cols2, &v2));
144:   }

146:   /* Test MatConvert() */
147:   PetscCall(MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C));

149:   /* See if MatMult Says both are same */
150:   for (i = 0; i < IMAX; i++) {
151:     PetscCall(VecSetRandom(xx, rand));
152:     PetscCall(MatMult(A, xx, s1));
153:     PetscCall(MatMult(C, xx, s2));
154:     PetscCall(VecNorm(s1, NORM_2, &s1norm));
155:     PetscCall(VecNorm(s2, NORM_2, &s2norm));
156:     rnorm = s2norm - s1norm;
157:     if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
158:   }
159:   PetscCall(MatDestroy(&C));

161:   /* Test MatTranspose() */
162:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At));
163:   PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &Bt));
164:   for (i = 0; i < IMAX; i++) {
165:     PetscCall(VecSetRandom(xx, rand));
166:     PetscCall(MatMult(At, xx, s1));
167:     PetscCall(MatMult(Bt, xx, s2));
168:     PetscCall(VecNorm(s1, NORM_2, &s1norm));
169:     PetscCall(VecNorm(s2, NORM_2, &s2norm));
170:     rnorm = s2norm - s1norm;
171:     if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs));
172:   }
173:   PetscCall(MatDestroy(&At));
174:   PetscCall(MatDestroy(&Bt));

176:   PetscCall(MatDestroy(&A));
177:   PetscCall(MatDestroy(&B));
178:   PetscCall(VecDestroy(&xx));
179:   PetscCall(VecDestroy(&yy));
180:   PetscCall(VecDestroy(&s1));
181:   PetscCall(VecDestroy(&s2));
182:   PetscCall(PetscRandomDestroy(&rand));
183:   PetscCall(PetscFinalize());
184:   return 0;
185: }

187: /*TEST

189:    build:
190:       requires: !complex

192:    test:
193:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
194:       nsize: 3
195:       args: -matload_block_size 1 -f ${DATAFILESPATH}/matrices/small

197:    test:
198:       suffix: 2
199:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
200:       nsize: 3
201:       args: -matload_block_size 2 -f ${DATAFILESPATH}/matrices/small
202:       output_file: output/ex53_1.out

204:    test:
205:       suffix: 3
206:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
207:       nsize: 3
208:       args: -matload_block_size 4 -f ${DATAFILESPATH}/matrices/small
209:       output_file: output/ex53_1.out

211:    test:
212:       TODO: Matrix row/column sizes are not compatible with block size
213:       suffix: 4
214:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
215:       nsize: 3
216:       args: -matload_block_size 5 -f ${DATAFILESPATH}/matrices/small
217:       output_file: output/ex53_1.out

219:    test:
220:       suffix: 5
221:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
222:       nsize: 3
223:       args: -matload_block_size 6 -f ${DATAFILESPATH}/matrices/small
224:       output_file: output/ex53_1.out

226:    test:
227:       TODO: Matrix row/column sizes are not compatible with block size
228:       suffix: 6
229:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
230:       nsize: 3
231:       args: -matload_block_size 7 -f ${DATAFILESPATH}/matrices/small
232:       output_file: output/ex53_1.out

234:    test:
235:       TODO: Matrix row/column sizes are not compatible with block size
236:       suffix: 7
237:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
238:       nsize: 3
239:       args: -matload_block_size 8 -f ${DATAFILESPATH}/matrices/small
240:       output_file: output/ex53_1.out

242:    test:
243:       suffix: 8
244:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
245:       nsize: 4
246:       args: -matload_block_size 3 -f ${DATAFILESPATH}/matrices/small
247:       output_file: output/ex53_1.out

249: TEST*/