Actual source code: ex53.c


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

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

 22:   PetscInitialize(&argc, &args, (char *)0, help);
 23:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 24:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 26:   /* Check out if MatLoad() works */
 27:   PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg);
 29:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
 30:   MatCreate(PETSC_COMM_WORLD, &A);
 31:   MatSetType(A, MATBAIJ);
 32:   MatLoad(A, fd);

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

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

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

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

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

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

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

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

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

123:     for (j = 0; j < 4; j++) {
124:       if (vals1[j] != vals2[j]) {
125:         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);
126:       }
127:     }
128:   }

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

137:     for (i = 0, j = 0; i < ncols1 && j < ncols2; j++) {
138:       while (cols2[j] != cols1[i]) i++;
140:     }

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

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

150:   /* See if MatMult Says both are same */
151:   for (i = 0; i < IMAX; i++) {
152:     VecSetRandom(xx, rand);
153:     MatMult(A, xx, s1);
154:     MatMult(C, xx, s2);
155:     VecNorm(s1, NORM_2, &s1norm);
156:     VecNorm(s2, NORM_2, &s2norm);
157:     rnorm = s2norm - s1norm;
158:     if (rnorm < -tol || rnorm > tol) 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);
159:   }
160:   MatDestroy(&C);

162:   /* Test MatTranspose() */
163:   MatTranspose(A, MAT_INITIAL_MATRIX, &At);
164:   MatTranspose(B, MAT_INITIAL_MATRIX, &Bt);
165:   for (i = 0; i < IMAX; i++) {
166:     VecSetRandom(xx, rand);
167:     MatMult(At, xx, s1);
168:     MatMult(Bt, xx, s2);
169:     VecNorm(s1, NORM_2, &s1norm);
170:     VecNorm(s2, NORM_2, &s2norm);
171:     rnorm = s2norm - s1norm;
172:     if (rnorm < -tol || rnorm > tol) 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);
173:   }
174:   MatDestroy(&At);
175:   MatDestroy(&Bt);

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

188: /*TEST

190:    build:
191:       requires: !complex

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

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

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

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

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

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

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

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

250: TEST*/