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