1 2 static char help[] = "Tests the various sequential routines in MatSBAIJ format. Same as ex74.c except diagonal entries of the matrices are zeros.\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Vec x, y, b, s1, s2; 9 Mat A; /* linear system matrix */ 10 Mat sA; /* symmetric part of the matrices */ 11 PetscInt n, mbs = 16, bs = 1, nz = 3, prob = 2, i, j, col[3], row, Ii, J, n1; 12 const PetscInt *ip_ptr; 13 PetscScalar neg_one = -1.0, value[3], alpha = 0.1; 14 PetscMPIInt size; 15 IS ip, isrow, iscol; 16 PetscRandom rdm; 17 PetscBool reorder = PETSC_FALSE; 18 MatInfo minfo1, minfo2; 19 PetscReal norm1, norm2, tol = 10 * PETSC_SMALL; 20 21 PetscFunctionBeginUser; 22 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 24 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL)); 26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL)); 27 28 n = mbs * bs; 29 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &A)); 30 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 31 PetscCall(MatCreateSeqSBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &sA)); 32 PetscCall(MatSetOption(sA, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 33 34 /* Test MatGetOwnershipRange() */ 35 PetscCall(MatGetOwnershipRange(A, &Ii, &J)); 36 PetscCall(MatGetOwnershipRange(sA, &i, &j)); 37 if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetOwnershipRange() in MatSBAIJ format\n")); 38 39 /* Assemble matrix */ 40 if (bs == 1) { 41 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL)); 42 if (prob == 1) { /* tridiagonal matrix */ 43 value[0] = -1.0; 44 value[1] = 2.0; 45 value[2] = -1.0; 46 for (i = 1; i < n - 1; i++) { 47 col[0] = i - 1; 48 col[1] = i; 49 col[2] = i + 1; 50 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 51 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES)); 52 } 53 i = n - 1; 54 col[0] = 0; 55 col[1] = n - 2; 56 col[2] = n - 1; 57 58 value[0] = 0.1; 59 value[1] = -1; 60 value[2] = 2; 61 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 62 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES)); 63 64 i = 0; 65 col[0] = 0; 66 col[1] = 1; 67 col[2] = n - 1; 68 69 value[0] = 2.0; 70 value[1] = -1.0; 71 value[2] = 0.1; 72 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 73 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES)); 74 } else if (prob == 2) { /* matrix for the five point stencil */ 75 n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001); 76 PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!"); 77 for (i = 0; i < n1; i++) { 78 for (j = 0; j < n1; j++) { 79 Ii = j + n1 * i; 80 if (i > 0) { 81 J = Ii - n1; 82 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 83 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 84 } 85 if (i < n1 - 1) { 86 J = Ii + n1; 87 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 88 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 89 } 90 if (j > 0) { 91 J = Ii - 1; 92 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 93 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 94 } 95 if (j < n1 - 1) { 96 J = Ii + 1; 97 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 98 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 99 } 100 } 101 } 102 } 103 } else { /* bs > 1 */ 104 #if defined(DIAGB) 105 for (block = 0; block < n / bs; block++) { 106 /* diagonal blocks */ 107 value[0] = -1.0; 108 value[1] = 4.0; 109 value[2] = -1.0; 110 for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) { 111 col[0] = i - 1; 112 col[1] = i; 113 col[2] = i + 1; 114 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 115 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES)); 116 } 117 i = bs - 1 + block * bs; 118 col[0] = bs - 2 + block * bs; 119 col[1] = bs - 1 + block * bs; 120 121 value[0] = -1.0; 122 value[1] = 4.0; 123 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); 124 PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES)); 125 126 i = 0 + block * bs; 127 col[0] = 0 + block * bs; 128 col[1] = 1 + block * bs; 129 130 value[0] = 4.0; 131 value[1] = -1.0; 132 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); 133 PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES)); 134 } 135 #endif 136 /* off-diagonal blocks */ 137 value[0] = -1.0; 138 for (i = 0; i < (n / bs - 1) * bs; i++) { 139 col[0] = i + bs; 140 PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES)); 141 PetscCall(MatSetValues(sA, 1, &i, 1, col, value, INSERT_VALUES)); 142 col[0] = i; 143 row = i + bs; 144 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES)); 145 PetscCall(MatSetValues(sA, 1, &row, 1, col, value, INSERT_VALUES)); 146 } 147 } 148 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 149 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 150 151 PetscCall(MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY)); 152 PetscCall(MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY)); 153 154 /* Test MatNorm() */ 155 PetscCall(MatNorm(A, NORM_FROBENIUS, &norm1)); 156 PetscCall(MatNorm(sA, NORM_FROBENIUS, &norm2)); 157 norm1 -= norm2; 158 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm(), fnorm1-fnorm2=%16.14e\n", (double)norm1)); 159 PetscCall(MatNorm(A, NORM_INFINITY, &norm1)); 160 PetscCall(MatNorm(sA, NORM_INFINITY, &norm2)); 161 norm1 -= norm2; 162 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n", (double)norm1)); 163 164 /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */ 165 PetscCall(MatGetInfo(A, MAT_LOCAL, &minfo1)); 166 PetscCall(MatGetInfo(sA, MAT_LOCAL, &minfo2)); 167 i = (int)(minfo1.nz_used - minfo2.nz_used); 168 j = (int)(minfo1.nz_allocated - minfo2.nz_allocated); 169 if (i < 0 || j < 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetInfo()\n")); 170 171 PetscCall(MatGetSize(A, &Ii, &J)); 172 PetscCall(MatGetSize(sA, &i, &j)); 173 if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetSize()\n")); 174 175 PetscCall(MatGetBlockSize(A, &Ii)); 176 PetscCall(MatGetBlockSize(sA, &i)); 177 if (i - Ii) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetBlockSize()\n")); 178 179 /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */ 180 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 181 PetscCall(PetscRandomSetFromOptions(rdm)); 182 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x)); 183 PetscCall(VecDuplicate(x, &s1)); 184 PetscCall(VecDuplicate(x, &s2)); 185 PetscCall(VecDuplicate(x, &y)); 186 PetscCall(VecDuplicate(x, &b)); 187 188 PetscCall(VecSetRandom(x, rdm)); 189 190 PetscCall(MatDiagonalScale(A, x, x)); 191 PetscCall(MatDiagonalScale(sA, x, x)); 192 193 PetscCall(MatGetDiagonal(A, s1)); 194 PetscCall(MatGetDiagonal(sA, s2)); 195 PetscCall(VecNorm(s1, NORM_1, &norm1)); 196 PetscCall(VecNorm(s2, NORM_1, &norm2)); 197 norm1 -= norm2; 198 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatGetDiagonal() \n")); 199 200 PetscCall(MatScale(A, alpha)); 201 PetscCall(MatScale(sA, alpha)); 202 203 /* Test MatMult(), MatMultAdd() */ 204 for (i = 0; i < 40; i++) { 205 PetscCall(VecSetRandom(x, rdm)); 206 PetscCall(MatMult(A, x, s1)); 207 PetscCall(MatMult(sA, x, s2)); 208 PetscCall(VecNorm(s1, NORM_1, &norm1)); 209 PetscCall(VecNorm(s2, NORM_1, &norm2)); 210 norm1 -= norm2; 211 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMult(), MatDiagonalScale() or MatScale()\n")); 212 } 213 214 for (i = 0; i < 40; i++) { 215 PetscCall(VecSetRandom(x, rdm)); 216 PetscCall(VecSetRandom(y, rdm)); 217 PetscCall(MatMultAdd(A, x, y, s1)); 218 PetscCall(MatMultAdd(sA, x, y, s2)); 219 PetscCall(VecNorm(s1, NORM_1, &norm1)); 220 PetscCall(VecNorm(s2, NORM_1, &norm2)); 221 norm1 -= norm2; 222 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n")); 223 } 224 225 /* Test MatReordering() */ 226 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &isrow, &iscol)); 227 ip = isrow; 228 229 if (reorder) { 230 IS nip; 231 PetscInt *nip_ptr; 232 PetscCall(PetscMalloc1(mbs, &nip_ptr)); 233 PetscCall(ISGetIndices(ip, &ip_ptr)); 234 PetscCall(PetscArraycpy(nip_ptr, ip_ptr, mbs)); 235 i = nip_ptr[1]; 236 nip_ptr[1] = nip_ptr[mbs - 2]; 237 nip_ptr[mbs - 2] = i; 238 i = nip_ptr[0]; 239 nip_ptr[0] = nip_ptr[mbs - 1]; 240 nip_ptr[mbs - 1] = i; 241 PetscCall(ISRestoreIndices(ip, &ip_ptr)); 242 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, mbs, nip_ptr, PETSC_COPY_VALUES, &nip)); 243 PetscCall(PetscFree(nip_ptr)); 244 245 PetscCall(MatReorderingSeqSBAIJ(sA, ip)); 246 PetscCall(ISDestroy(&nip)); 247 } 248 249 PetscCall(ISDestroy(&iscol)); 250 PetscCall(ISDestroy(&isrow)); 251 PetscCall(MatDestroy(&A)); 252 PetscCall(MatDestroy(&sA)); 253 PetscCall(VecDestroy(&x)); 254 PetscCall(VecDestroy(&y)); 255 PetscCall(VecDestroy(&s1)); 256 PetscCall(VecDestroy(&s2)); 257 PetscCall(VecDestroy(&b)); 258 PetscCall(PetscRandomDestroy(&rdm)); 259 260 PetscCall(PetscFinalize()); 261 return 0; 262 } 263 264 /*TEST 265 266 test: 267 args: -bs {{1 2 3 4 5 6 7 8}} 268 269 TEST*/ 270