Actual source code: ex48.c

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

  3: #include <petscmat.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

155:   PetscCall(MatFactorInfoInitialize(&info));

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

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

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

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

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

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

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

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

243: /*TEST

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

248: TEST*/