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