1 static char help[] = "Tests MatGetRowMax(), MatGetRowMin(), MatGetRowMaxAbs()\n"; 2 3 #include <petscmat.h> 4 5 #define M 5 6 #define N 6 7 8 int main(int argc, char **args) 9 { 10 Mat A; 11 Vec min, max, maxabs, e; 12 PetscInt m, n, j, imin[M], imax[M], imaxabs[M], indices[N], row, testcase = 0; 13 PetscScalar values[N]; 14 PetscMPIInt size, rank; 15 PetscReal enorm; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 20 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-testcase", &testcase, NULL)); 22 23 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 24 if (testcase == 1) { /* proc[0] holds entire A and other processes have no entry */ 25 if (rank == 0) { 26 PetscCall(MatSetSizes(A, M, N, PETSC_DECIDE, PETSC_DECIDE)); 27 } else { 28 PetscCall(MatSetSizes(A, 0, 0, PETSC_DECIDE, PETSC_DECIDE)); 29 } 30 testcase = 0; 31 } else { 32 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N)); 33 } 34 PetscCall(MatSetFromOptions(A)); 35 PetscCall(MatSetUp(A)); 36 37 if (rank == 0) { /* proc[0] sets matrix A */ 38 for (j = 0; j < N; j++) indices[j] = j; 39 switch (testcase) { 40 case 1: /* see testcast 0 */ 41 break; 42 case 2: 43 row = 0; 44 values[0] = -2.0; 45 values[1] = -2.0; 46 values[2] = -2.0; 47 values[3] = -4.0; 48 values[4] = 1.0; 49 values[5] = 1.0; 50 PetscCall(MatSetValues(A, 1, &row, N, indices, values, INSERT_VALUES)); 51 row = 2; 52 indices[0] = 0; 53 indices[1] = 3; 54 indices[2] = 5; 55 values[0] = -2.0; 56 values[1] = -2.0; 57 values[2] = -2.0; 58 PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES)); 59 row = 3; 60 indices[0] = 0; 61 indices[1] = 1; 62 indices[2] = 4; 63 values[0] = -2.0; 64 values[1] = -2.0; 65 values[2] = -2.0; 66 PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES)); 67 row = 4; 68 indices[0] = 0; 69 indices[1] = 1; 70 indices[2] = 2; 71 values[0] = -2.0; 72 values[1] = -2.0; 73 values[2] = -2.0; 74 PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES)); 75 break; 76 case 3: 77 row = 0; 78 values[0] = -2.0; 79 values[1] = -2.0; 80 values[2] = -2.0; 81 PetscCall(MatSetValues(A, 1, &row, 3, indices + 1, values, INSERT_VALUES)); 82 row = 1; 83 values[0] = -2.0; 84 values[1] = -2.0; 85 values[2] = -2.0; 86 PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES)); 87 row = 2; 88 values[0] = -2.0; 89 values[1] = -2.0; 90 values[2] = -2.0; 91 PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES)); 92 row = 3; 93 values[0] = -2.0; 94 values[1] = -2.0; 95 values[2] = -2.0; 96 values[3] = -1.0; 97 PetscCall(MatSetValues(A, 1, &row, 4, indices, values, INSERT_VALUES)); 98 row = 4; 99 values[0] = -2.0; 100 values[1] = -2.0; 101 values[2] = -2.0; 102 values[3] = -1.0; 103 PetscCall(MatSetValues(A, 1, &row, 4, indices, values, INSERT_VALUES)); 104 break; 105 106 default: 107 row = 0; 108 values[0] = -1.0; 109 values[1] = 0.0; 110 values[2] = 1.0; 111 values[3] = 3.0; 112 values[4] = 4.0; 113 values[5] = -5.0; 114 PetscCall(MatSetValues(A, 1, &row, N, indices, values, INSERT_VALUES)); 115 row = 1; 116 PetscCall(MatSetValues(A, 1, &row, 3, indices, values, INSERT_VALUES)); 117 row = 3; 118 PetscCall(MatSetValues(A, 1, &row, 1, indices + 4, values + 4, INSERT_VALUES)); 119 row = 4; 120 PetscCall(MatSetValues(A, 1, &row, 2, indices + 4, values + 4, INSERT_VALUES)); 121 } 122 } 123 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 124 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 125 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 126 127 PetscCall(MatGetLocalSize(A, &m, &n)); 128 PetscCall(VecCreate(PETSC_COMM_WORLD, &min)); 129 PetscCall(VecSetSizes(min, m, PETSC_DECIDE)); 130 PetscCall(VecSetFromOptions(min)); 131 PetscCall(VecDuplicate(min, &max)); 132 PetscCall(VecDuplicate(min, &maxabs)); 133 PetscCall(VecDuplicate(min, &e)); 134 135 /* MatGetRowMax() */ 136 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMax\n")); 137 PetscCall(MatGetRowMax(A, max, NULL)); 138 PetscCall(MatGetRowMax(A, max, imax)); 139 PetscCall(VecView(max, PETSC_VIEWER_STDOUT_WORLD)); 140 PetscCall(VecGetLocalSize(max, &n)); 141 PetscCall(PetscIntView(n, imax, PETSC_VIEWER_STDOUT_WORLD)); 142 143 /* MatGetRowMin() */ 144 PetscCall(MatScale(A, -1.0)); 145 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 146 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMin\n")); 147 PetscCall(MatGetRowMin(A, min, NULL)); 148 PetscCall(MatGetRowMin(A, min, imin)); 149 150 PetscCall(VecWAXPY(e, 1.0, max, min)); /* e = max + min */ 151 PetscCall(VecNorm(e, NORM_INFINITY, &enorm)); 152 PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "max+min > PETSC_MACHINE_EPSILON "); 153 for (j = 0; j < n; j++) PetscCheck(imin[j] == imax[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT, j, imin[j], imax[j]); 154 155 /* MatGetRowMaxAbs() */ 156 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMaxAbs\n")); 157 PetscCall(MatGetRowMaxAbs(A, maxabs, NULL)); 158 PetscCall(MatGetRowMaxAbs(A, maxabs, imaxabs)); 159 PetscCall(VecView(maxabs, PETSC_VIEWER_STDOUT_WORLD)); 160 PetscCall(PetscIntView(n, imaxabs, PETSC_VIEWER_STDOUT_WORLD)); 161 162 /* MatGetRowMinAbs() */ 163 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMinAbs\n")); 164 PetscCall(MatGetRowMinAbs(A, min, NULL)); 165 PetscCall(MatGetRowMinAbs(A, min, imin)); 166 PetscCall(VecView(min, PETSC_VIEWER_STDOUT_WORLD)); 167 PetscCall(PetscIntView(n, imin, PETSC_VIEWER_STDOUT_WORLD)); 168 169 if (size == 1) { 170 /* Test MatGetRowMax, MatGetRowMin and MatGetRowMaxAbs for SeqDense and MPIBAIJ matrix */ 171 Mat Adense; 172 Vec max_d, maxabs_d; 173 PetscCall(VecDuplicate(min, &max_d)); 174 PetscCall(VecDuplicate(min, &maxabs_d)); 175 176 PetscCall(MatScale(A, -1.0)); 177 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense)); 178 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatGetRowMax for seqdense matrix\n")); 179 PetscCall(MatGetRowMax(Adense, max_d, imax)); 180 181 PetscCall(VecWAXPY(e, -1.0, max, max_d)); /* e = -max + max_d */ 182 PetscCall(VecNorm(e, NORM_INFINITY, &enorm)); 183 PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "norm(-max + max_d) %g > PETSC_MACHINE_EPSILON", (double)enorm); 184 185 PetscCall(MatScale(Adense, -1.0)); 186 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatGetRowMin for seqdense matrix\n")); 187 PetscCall(MatGetRowMin(Adense, min, imin)); 188 189 PetscCall(VecWAXPY(e, 1.0, max, min)); /* e = max + min */ 190 PetscCall(VecNorm(e, NORM_INFINITY, &enorm)); 191 PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "max+min > PETSC_MACHINE_EPSILON "); 192 for (j = 0; j < n; j++) PetscCheck(imin[j] == imax[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT " for seqdense matrix", j, imin[j], imax[j]); 193 194 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatGetRowMaxAbs for seqdense matrix\n")); 195 PetscCall(MatGetRowMaxAbs(Adense, maxabs_d, imaxabs)); 196 PetscCall(VecWAXPY(e, -1.0, maxabs, maxabs_d)); /* e = -maxabs + maxabs_d */ 197 PetscCall(VecNorm(e, NORM_INFINITY, &enorm)); 198 PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON", (double)enorm); 199 200 PetscCall(MatDestroy(&Adense)); 201 PetscCall(VecDestroy(&max_d)); 202 PetscCall(VecDestroy(&maxabs_d)); 203 } 204 205 { /* BAIJ matrix */ 206 Mat B; 207 Vec maxabsB, maxabsB2; 208 PetscInt bs = 2, *imaxabsB, *imaxabsB2, rstart, rend, cstart, cend, ncols, col, Brows[2], Bcols[2]; 209 const PetscInt *cols; 210 const PetscScalar *vals, *vals2; 211 PetscScalar Bvals[4]; 212 213 PetscCall(PetscMalloc2(M, &imaxabsB, bs * M, &imaxabsB2)); 214 215 /* bs = 1 */ 216 PetscCall(MatConvert(A, MATMPIBAIJ, MAT_INITIAL_MATRIX, &B)); 217 PetscCall(VecDuplicate(min, &maxabsB)); 218 PetscCall(MatGetRowMaxAbs(B, maxabsB, NULL)); 219 PetscCall(MatGetRowMaxAbs(B, maxabsB, imaxabsB)); 220 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MatGetRowMaxAbs for BAIJ matrix\n")); 221 PetscCall(VecWAXPY(e, -1.0, maxabs, maxabsB)); /* e = -maxabs + maxabsB */ 222 PetscCall(VecNorm(e, NORM_INFINITY, &enorm)); 223 PetscCheck(enorm <= PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON", (double)enorm); 224 225 for (j = 0; j < n; j++) PetscCheck(imaxabs[j] == imaxabsB[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "imaxabs[%" PetscInt_FMT "] %" PetscInt_FMT " != imaxabsB %" PetscInt_FMT, j, imin[j], imax[j]); 226 PetscCall(MatDestroy(&B)); 227 228 /* Test bs = 2: Create B with bs*bs block structure of A */ 229 PetscCall(VecCreate(PETSC_COMM_WORLD, &maxabsB2)); 230 PetscCall(VecSetSizes(maxabsB2, bs * m, PETSC_DECIDE)); 231 PetscCall(VecSetFromOptions(maxabsB2)); 232 233 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 234 PetscCall(MatGetOwnershipRangeColumn(A, &cstart, &cend)); 235 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 236 PetscCall(MatSetSizes(B, bs * (rend - rstart), bs * (cend - cstart), PETSC_DECIDE, PETSC_DECIDE)); 237 PetscCall(MatSetFromOptions(B)); 238 PetscCall(MatSetUp(B)); 239 240 for (row = rstart; row < rend; row++) { 241 PetscCall(MatGetRow(A, row, &ncols, &cols, &vals)); 242 for (col = 0; col < ncols; col++) { 243 for (j = 0; j < bs; j++) { 244 Brows[j] = bs * row + j; 245 Bcols[j] = bs * cols[col] + j; 246 } 247 for (j = 0; j < bs * bs; j++) Bvals[j] = vals[col]; 248 PetscCall(MatSetValues(B, bs, Brows, bs, Bcols, Bvals, INSERT_VALUES)); 249 } 250 PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals)); 251 } 252 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 253 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 254 255 PetscCall(MatGetRowMaxAbs(B, maxabsB2, imaxabsB2)); 256 257 /* Check maxabsB2 and imaxabsB2 */ 258 PetscCall(VecGetArrayRead(maxabsB, &vals)); 259 PetscCall(VecGetArrayRead(maxabsB2, &vals2)); 260 for (row = 0; row < m; row++) PetscCheck(PetscAbsScalar(vals[row] - vals2[bs * row]) <= PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_PLIB, "row %" PetscInt_FMT " maxabsB != maxabsB2", row); 261 PetscCall(VecRestoreArrayRead(maxabsB, &vals)); 262 PetscCall(VecRestoreArrayRead(maxabsB2, &vals2)); 263 264 for (col = 0; col < n; col++) PetscCheck(imaxabsB[col] == imaxabsB2[bs * col] / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "col %" PetscInt_FMT " imaxabsB != imaxabsB2", col); 265 PetscCall(VecDestroy(&maxabsB)); 266 PetscCall(MatDestroy(&B)); 267 PetscCall(VecDestroy(&maxabsB2)); 268 PetscCall(PetscFree2(imaxabsB, imaxabsB2)); 269 } 270 271 PetscCall(VecDestroy(&min)); 272 PetscCall(VecDestroy(&max)); 273 PetscCall(VecDestroy(&maxabs)); 274 PetscCall(VecDestroy(&e)); 275 PetscCall(MatDestroy(&A)); 276 PetscCall(PetscFinalize()); 277 return 0; 278 } 279 280 /*TEST 281 282 test: 283 output_file: output/ex114.out 284 285 test: 286 suffix: 2 287 args: -testcase 1 288 output_file: output/ex114.out 289 290 test: 291 suffix: 3 292 args: -testcase 2 293 output_file: output/ex114_3.out 294 295 test: 296 suffix: 4 297 args: -testcase 3 298 output_file: output/ex114_4.out 299 300 test: 301 suffix: 5 302 nsize: 3 303 args: -testcase 0 304 output_file: output/ex114_5.out 305 306 test: 307 suffix: 6 308 nsize: 3 309 args: -testcase 1 310 output_file: output/ex114_6.out 311 312 test: 313 suffix: 7 314 nsize: 3 315 args: -testcase 2 316 output_file: output/ex114_7.out 317 318 test: 319 suffix: 8 320 nsize: 3 321 args: -testcase 3 322 output_file: output/ex114_8.out 323 324 TEST*/ 325