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