1 static char help[] = "Tests various routines in MatSeqBAIJ format.\n"; 2 3 #include <petscmat.h> 4 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; 16 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)); 33 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 } 43 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 } 61 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)); 66 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)); 80 81 /* MatShift() */ 82 rval = 10 * s1norm; 83 PetscCall(MatShift(A, rval)); 84 PetscCall(MatShift(B, rval)); 85 86 /* Test MatTranspose() */ 87 PetscCall(MatTranspose(A, MAT_INPLACE_MATRIX, &A)); 88 PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B)); 89 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 } 105 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 } 117 118 /* Test MatMult() */ 119 PetscCall(MatMultEqual(A, B, 10, &flg)); 120 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMult()\n")); 121 122 /* Test MatMultAdd() */ 123 PetscCall(MatMultAddEqual(A, B, 10, &flg)); 124 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultAdd()\n")); 125 126 /* Test MatMultTranspose() */ 127 PetscCall(MatMultTransposeEqual(A, B, 10, &flg)); 128 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTranspose()\n")); 129 130 /* Test MatMultTransposeAdd() */ 131 PetscCall(MatMultTransposeAddEqual(A, B, 10, &flg)); 132 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTransposeAdd()\n")); 133 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")); 145 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)); 154 155 PetscCall(MatFactorInfoInitialize(&info)); 156 157 info.fill = 2.0; 158 info.dtcol = 0.0; 159 info.zeropivot = 1.e-14; 160 info.pivotinblocks = 1.0; 161 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 } 175 176 PetscCall(MatLUFactor(B, is1, is2, &info)); 177 PetscCall(MatLUFactor(A, is1, is2, &info)); 178 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 } 190 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 } 203 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 } 214 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 } 227 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 } 242 243 /*TEST 244 245 test: 246 args: -mat_block_size {{1 2 3 4 5 6 7 8}} 247 248 TEST*/ 249