Actual source code: ex48.c


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

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Mat           A, B, C, D, Fact;
  9:   Vec           xx, s1, s2, yy;
 10:   PetscInt      m = 45, rows[2], cols[2], bs = 1, i, row, col, *idx, M;
 11:   PetscScalar   rval, vals1[4], vals2[4];
 12:   PetscRandom   rdm;
 13:   IS            is1, is2;
 14:   PetscReal     s1norm, s2norm, rnorm, tol = 1.e-4;
 15:   PetscBool     flg;
 16:   MatFactorInfo info;

 19:   PetscInitialize(&argc, &args, (char *)0, help);
 20:   /* Test MatSetValues() and MatGetValues() */
 21:   PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL);
 22:   PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL);
 23:   M = m * bs;
 24:   MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A);
 25:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 26:   MatCreateSeqAIJ(PETSC_COMM_SELF, M, M, 15, NULL, &B);
 27:   MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 28:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
 29:   PetscRandomSetFromOptions(rdm);
 30:   VecCreateSeq(PETSC_COMM_SELF, M, &xx);
 31:   VecDuplicate(xx, &s1);
 32:   VecDuplicate(xx, &s2);
 33:   VecDuplicate(xx, &yy);

 35:   /* For each row add at least 15 elements */
 36:   for (row = 0; row < M; row++) {
 37:     for (i = 0; i < 25 * bs; i++) {
 38:       PetscRandomGetValue(rdm, &rval);
 39:       col = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
 40:       MatSetValues(A, 1, &row, 1, &col, &rval, INSERT_VALUES);
 41:       MatSetValues(B, 1, &row, 1, &col, &rval, INSERT_VALUES);
 42:     }
 43:   }

 45:   /* Now set blocks of values */
 46:   for (i = 0; i < 20 * bs; i++) {
 47:     PetscRandomGetValue(rdm, &rval);
 48:     cols[0]  = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
 49:     vals1[0] = rval;
 50:     PetscRandomGetValue(rdm, &rval);
 51:     cols[1]  = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
 52:     vals1[1] = rval;
 53:     PetscRandomGetValue(rdm, &rval);
 54:     rows[0]  = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
 55:     vals1[2] = rval;
 56:     PetscRandomGetValue(rdm, &rval);
 57:     rows[1]  = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
 58:     vals1[3] = rval;
 59:     MatSetValues(A, 2, rows, 2, cols, vals1, INSERT_VALUES);
 60:     MatSetValues(B, 2, rows, 2, cols, vals1, INSERT_VALUES);
 61:   }

 63:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 64:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 65:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

 68:   /* Test MatNorm() */
 69:   MatNorm(A, NORM_FROBENIUS, &s1norm);
 70:   MatNorm(B, NORM_FROBENIUS, &s2norm);
 71:   rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
 72:   if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
 73:   MatNorm(A, NORM_INFINITY, &s1norm);
 74:   MatNorm(B, NORM_INFINITY, &s2norm);
 75:   rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
 76:   if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
 77:   MatNorm(A, NORM_1, &s1norm);
 78:   MatNorm(B, NORM_1, &s2norm);
 79:   rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
 80:   if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);

 82:   /* MatShift() */
 83:   rval = 10 * s1norm;
 84:   MatShift(A, rval);
 85:   MatShift(B, rval);

 87:   /* Test MatTranspose() */
 88:   MatTranspose(A, MAT_INPLACE_MATRIX, &A);
 89:   MatTranspose(B, MAT_INPLACE_MATRIX, &B);

 91:   /* Now do MatGetValues()  */
 92:   for (i = 0; i < 30; i++) {
 93:     PetscRandomGetValue(rdm, &rval);
 94:     cols[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
 95:     PetscRandomGetValue(rdm, &rval);
 96:     cols[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
 97:     PetscRandomGetValue(rdm, &rval);
 98:     rows[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
 99:     PetscRandomGetValue(rdm, &rval);
100:     rows[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
101:     MatGetValues(A, 2, rows, 2, cols, vals1);
102:     MatGetValues(B, 2, rows, 2, cols, vals2);
103:     PetscArraycmp(vals1, vals2, 4, &flg);
104:     if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetValues bs = %" PetscInt_FMT "\n", bs);
105:   }

107:   /* Test MatMult(), MatMultAdd() */
108:   for (i = 0; i < 40; i++) {
109:     VecSetRandom(xx, rdm);
110:     VecSet(s2, 0.0);
111:     MatMult(A, xx, s1);
112:     MatMultAdd(A, xx, s2, s2);
113:     VecNorm(s1, NORM_2, &s1norm);
114:     VecNorm(s2, NORM_2, &s2norm);
115:     rnorm = s2norm - s1norm;
116:     if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
117:   }

119:   /* Test MatMult() */
120:   MatMultEqual(A, B, 10, &flg);
121:   if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMult()\n");

123:   /* Test MatMultAdd() */
124:   MatMultAddEqual(A, B, 10, &flg);
125:   if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMultAdd()\n");

127:   /* Test MatMultTranspose() */
128:   MatMultTransposeEqual(A, B, 10, &flg);
129:   if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTranspose()\n");

131:   /* Test MatMultTransposeAdd() */
132:   MatMultTransposeAddEqual(A, B, 10, &flg);
133:   if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTransposeAdd()\n");

135:   /* Test MatMatMult() */
136:   MatCreateSeqDense(PETSC_COMM_SELF, M, 40, NULL, &C);
137:   MatSetRandom(C, rdm);
138:   MatMatMult(A, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &D);
139:   MatMatMultEqual(A, C, D, 40, &flg);
140:   if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult()\n");
141:   MatDestroy(&D);
142:   MatCreateSeqDense(PETSC_COMM_SELF, M, 40, NULL, &D);
143:   MatMatMult(A, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &D);
144:   MatMatMultEqual(A, C, D, 40, &flg);
145:   if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult()\n");

147:   /* Do LUFactor() on both the matrices */
148:   PetscMalloc1(M, &idx);
149:   for (i = 0; i < M; i++) idx[i] = i;
150:   ISCreateGeneral(PETSC_COMM_SELF, M, idx, PETSC_COPY_VALUES, &is1);
151:   ISCreateGeneral(PETSC_COMM_SELF, M, idx, PETSC_COPY_VALUES, &is2);
152:   PetscFree(idx);
153:   ISSetPermutation(is1);
154:   ISSetPermutation(is2);

156:   MatFactorInfoInitialize(&info);

158:   info.fill          = 2.0;
159:   info.dtcol         = 0.0;
160:   info.zeropivot     = 1.e-14;
161:   info.pivotinblocks = 1.0;

163:   if (bs < 4) {
164:     MatGetFactor(A, "petsc", MAT_FACTOR_LU, &Fact);
165:     MatLUFactorSymbolic(Fact, A, is1, is2, &info);
166:     MatLUFactorNumeric(Fact, A, &info);
167:     VecSetRandom(yy, rdm);
168:     MatForwardSolve(Fact, yy, xx);
169:     MatBackwardSolve(Fact, xx, s1);
170:     MatDestroy(&Fact);
171:     VecScale(s1, -1.0);
172:     MatMultAdd(A, s1, yy, yy);
173:     VecNorm(yy, NORM_2, &rnorm);
174:     if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %" PetscInt_FMT "\n", (double)rnorm, bs);
175:   }

177:   MatLUFactor(B, is1, is2, &info);
178:   MatLUFactor(A, is1, is2, &info);

180:   /* Test MatSolveAdd() */
181:   for (i = 0; i < 10; i++) {
182:     VecSetRandom(xx, rdm);
183:     VecSetRandom(yy, rdm);
184:     MatSolveAdd(B, xx, yy, s2);
185:     MatSolveAdd(A, xx, yy, s1);
186:     VecNorm(s1, NORM_2, &s1norm);
187:     VecNorm(s2, NORM_2, &s2norm);
188:     rnorm = s2norm - s1norm;
189:     if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
190:   }

192:   /* Test MatSolveAdd() when x = A'b +x */
193:   for (i = 0; i < 10; i++) {
194:     VecSetRandom(xx, rdm);
195:     VecSetRandom(s1, rdm);
196:     VecCopy(s2, s1);
197:     MatSolveAdd(B, xx, s2, s2);
198:     MatSolveAdd(A, xx, s1, s1);
199:     VecNorm(s1, NORM_2, &s1norm);
200:     VecNorm(s2, NORM_2, &s2norm);
201:     rnorm = s2norm - s1norm;
202:     if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
203:   }

205:   /* Test MatSolve() */
206:   for (i = 0; i < 10; i++) {
207:     VecSetRandom(xx, rdm);
208:     MatSolve(B, xx, s2);
209:     MatSolve(A, xx, s1);
210:     VecNorm(s1, NORM_2, &s1norm);
211:     VecNorm(s2, NORM_2, &s2norm);
212:     rnorm = s2norm - s1norm;
213:     if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
214:   }

216:   /* Test MatSolveTranspose() */
217:   if (bs < 8) {
218:     for (i = 0; i < 10; i++) {
219:       VecSetRandom(xx, rdm);
220:       MatSolveTranspose(B, xx, s2);
221:       MatSolveTranspose(A, xx, s1);
222:       VecNorm(s1, NORM_2, &s1norm);
223:       VecNorm(s2, NORM_2, &s2norm);
224:       rnorm = s2norm - s1norm;
225:       if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
226:     }
227:   }

229:   MatDestroy(&A);
230:   MatDestroy(&B);
231:   MatDestroy(&C);
232:   MatDestroy(&D);
233:   VecDestroy(&xx);
234:   VecDestroy(&s1);
235:   VecDestroy(&s2);
236:   VecDestroy(&yy);
237:   ISDestroy(&is1);
238:   ISDestroy(&is2);
239:   PetscRandomDestroy(&rdm);
240:   PetscFinalize();
241:   return 0;
242: }

244: /*TEST

246:    test:
247:       args: -mat_block_size {{1 2 3 4 5 6 7 8}}

249: TEST*/