1 2 static char help[] = "Tests various routines in MatSeqBAIJ format.\n"; 3 4 #include <petscmat.h> 5 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; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 20 /* Test MatSetValues() and MatGetValues() */ 21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL)); 22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL)); 23 M = m * bs; 24 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A)); 25 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 26 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, M, M, 15, NULL, &B)); 27 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 28 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 29 PetscCall(PetscRandomSetFromOptions(rdm)); 30 PetscCall(VecCreateSeq(PETSC_COMM_SELF, M, &xx)); 31 PetscCall(VecDuplicate(xx, &s1)); 32 PetscCall(VecDuplicate(xx, &s2)); 33 PetscCall(VecDuplicate(xx, &yy)); 34 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 PetscCall(PetscRandomGetValue(rdm, &rval)); 39 col = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 40 PetscCall(MatSetValues(A, 1, &row, 1, &col, &rval, INSERT_VALUES)); 41 PetscCall(MatSetValues(B, 1, &row, 1, &col, &rval, INSERT_VALUES)); 42 } 43 } 44 45 /* Now set blocks of values */ 46 for (i = 0; i < 20 * bs; i++) { 47 PetscCall(PetscRandomGetValue(rdm, &rval)); 48 cols[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 49 vals1[0] = rval; 50 PetscCall(PetscRandomGetValue(rdm, &rval)); 51 cols[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 52 vals1[1] = rval; 53 PetscCall(PetscRandomGetValue(rdm, &rval)); 54 rows[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 55 vals1[2] = rval; 56 PetscCall(PetscRandomGetValue(rdm, &rval)); 57 rows[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 58 vals1[3] = rval; 59 PetscCall(MatSetValues(A, 2, rows, 2, cols, vals1, INSERT_VALUES)); 60 PetscCall(MatSetValues(B, 2, rows, 2, cols, vals1, INSERT_VALUES)); 61 } 62 63 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 64 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 65 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 66 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 67 68 /* Test MatNorm() */ 69 PetscCall(MatNorm(A, NORM_FROBENIUS, &s1norm)); 70 PetscCall(MatNorm(B, NORM_FROBENIUS, &s2norm)); 71 rnorm = PetscAbsReal(s2norm - s1norm) / s2norm; 72 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)); 73 PetscCall(MatNorm(A, NORM_INFINITY, &s1norm)); 74 PetscCall(MatNorm(B, NORM_INFINITY, &s2norm)); 75 rnorm = PetscAbsReal(s2norm - s1norm) / s2norm; 76 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)); 77 PetscCall(MatNorm(A, NORM_1, &s1norm)); 78 PetscCall(MatNorm(B, NORM_1, &s2norm)); 79 rnorm = PetscAbsReal(s2norm - s1norm) / s2norm; 80 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 82 /* MatShift() */ 83 rval = 10 * s1norm; 84 PetscCall(MatShift(A, rval)); 85 PetscCall(MatShift(B, rval)); 86 87 /* Test MatTranspose() */ 88 PetscCall(MatTranspose(A, MAT_INPLACE_MATRIX, &A)); 89 PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B)); 90 91 /* Now do MatGetValues() */ 92 for (i = 0; i < 30; i++) { 93 PetscCall(PetscRandomGetValue(rdm, &rval)); 94 cols[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 95 PetscCall(PetscRandomGetValue(rdm, &rval)); 96 cols[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 97 PetscCall(PetscRandomGetValue(rdm, &rval)); 98 rows[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 99 PetscCall(PetscRandomGetValue(rdm, &rval)); 100 rows[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M)); 101 PetscCall(MatGetValues(A, 2, rows, 2, cols, vals1)); 102 PetscCall(MatGetValues(B, 2, rows, 2, cols, vals2)); 103 PetscCall(PetscArraycmp(vals1, vals2, 4, &flg)); 104 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetValues bs = %" PetscInt_FMT "\n", bs)); 105 } 106 107 /* Test MatMult(), MatMultAdd() */ 108 for (i = 0; i < 40; i++) { 109 PetscCall(VecSetRandom(xx, rdm)); 110 PetscCall(VecSet(s2, 0.0)); 111 PetscCall(MatMult(A, xx, s1)); 112 PetscCall(MatMultAdd(A, xx, s2, s2)); 113 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 114 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 115 rnorm = s2norm - s1norm; 116 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)); 117 } 118 119 /* Test MatMult() */ 120 PetscCall(MatMultEqual(A, B, 10, &flg)); 121 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMult()\n")); 122 123 /* Test MatMultAdd() */ 124 PetscCall(MatMultAddEqual(A, B, 10, &flg)); 125 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultAdd()\n")); 126 127 /* Test MatMultTranspose() */ 128 PetscCall(MatMultTransposeEqual(A, B, 10, &flg)); 129 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTranspose()\n")); 130 131 /* Test MatMultTransposeAdd() */ 132 PetscCall(MatMultTransposeAddEqual(A, B, 10, &flg)); 133 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTransposeAdd()\n")); 134 135 /* Test MatMatMult() */ 136 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, M, 40, NULL, &C)); 137 PetscCall(MatSetRandom(C, rdm)); 138 PetscCall(MatMatMult(A, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &D)); 139 PetscCall(MatMatMultEqual(A, C, D, 40, &flg)); 140 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult()\n")); 141 PetscCall(MatDestroy(&D)); 142 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, M, 40, NULL, &D)); 143 PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &D)); 144 PetscCall(MatMatMultEqual(A, C, D, 40, &flg)); 145 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult()\n")); 146 147 /* Do LUFactor() on both the matrices */ 148 PetscCall(PetscMalloc1(M, &idx)); 149 for (i = 0; i < M; i++) idx[i] = i; 150 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, M, idx, PETSC_COPY_VALUES, &is1)); 151 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, M, idx, PETSC_COPY_VALUES, &is2)); 152 PetscCall(PetscFree(idx)); 153 PetscCall(ISSetPermutation(is1)); 154 PetscCall(ISSetPermutation(is2)); 155 156 PetscCall(MatFactorInfoInitialize(&info)); 157 158 info.fill = 2.0; 159 info.dtcol = 0.0; 160 info.zeropivot = 1.e-14; 161 info.pivotinblocks = 1.0; 162 163 if (bs < 4) { 164 PetscCall(MatGetFactor(A, "petsc", MAT_FACTOR_LU, &Fact)); 165 PetscCall(MatLUFactorSymbolic(Fact, A, is1, is2, &info)); 166 PetscCall(MatLUFactorNumeric(Fact, A, &info)); 167 PetscCall(VecSetRandom(yy, rdm)); 168 PetscCall(MatForwardSolve(Fact, yy, xx)); 169 PetscCall(MatBackwardSolve(Fact, xx, s1)); 170 PetscCall(MatDestroy(&Fact)); 171 PetscCall(VecScale(s1, -1.0)); 172 PetscCall(MatMultAdd(A, s1, yy, yy)); 173 PetscCall(VecNorm(yy, NORM_2, &rnorm)); 174 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %" PetscInt_FMT "\n", (double)rnorm, bs)); 175 } 176 177 PetscCall(MatLUFactor(B, is1, is2, &info)); 178 PetscCall(MatLUFactor(A, is1, is2, &info)); 179 180 /* Test MatSolveAdd() */ 181 for (i = 0; i < 10; i++) { 182 PetscCall(VecSetRandom(xx, rdm)); 183 PetscCall(VecSetRandom(yy, rdm)); 184 PetscCall(MatSolveAdd(B, xx, yy, s2)); 185 PetscCall(MatSolveAdd(A, xx, yy, s1)); 186 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 187 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 188 rnorm = s2norm - s1norm; 189 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)); 190 } 191 192 /* Test MatSolveAdd() when x = A'b +x */ 193 for (i = 0; i < 10; i++) { 194 PetscCall(VecSetRandom(xx, rdm)); 195 PetscCall(VecSetRandom(s1, rdm)); 196 PetscCall(VecCopy(s2, s1)); 197 PetscCall(MatSolveAdd(B, xx, s2, s2)); 198 PetscCall(MatSolveAdd(A, xx, s1, s1)); 199 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 200 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 201 rnorm = s2norm - s1norm; 202 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)); 203 } 204 205 /* Test MatSolve() */ 206 for (i = 0; i < 10; i++) { 207 PetscCall(VecSetRandom(xx, rdm)); 208 PetscCall(MatSolve(B, xx, s2)); 209 PetscCall(MatSolve(A, xx, s1)); 210 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 211 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 212 rnorm = s2norm - s1norm; 213 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)); 214 } 215 216 /* Test MatSolveTranspose() */ 217 if (bs < 8) { 218 for (i = 0; i < 10; i++) { 219 PetscCall(VecSetRandom(xx, rdm)); 220 PetscCall(MatSolveTranspose(B, xx, s2)); 221 PetscCall(MatSolveTranspose(A, xx, s1)); 222 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 223 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 224 rnorm = s2norm - s1norm; 225 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)); 226 } 227 } 228 229 PetscCall(MatDestroy(&A)); 230 PetscCall(MatDestroy(&B)); 231 PetscCall(MatDestroy(&C)); 232 PetscCall(MatDestroy(&D)); 233 PetscCall(VecDestroy(&xx)); 234 PetscCall(VecDestroy(&s1)); 235 PetscCall(VecDestroy(&s2)); 236 PetscCall(VecDestroy(&yy)); 237 PetscCall(ISDestroy(&is1)); 238 PetscCall(ISDestroy(&is2)); 239 PetscCall(PetscRandomDestroy(&rdm)); 240 PetscCall(PetscFinalize()); 241 return 0; 242 } 243 244 /*TEST 245 246 test: 247 args: -mat_block_size {{1 2 3 4 5 6 7 8}} 248 249 TEST*/ 250