1 2 static char help[] = "Tests various routines in MatMPIBAIJ format.\n"; 3 4 #include <petscmat.h> 5 #define IMAX 15 6 int main(int argc, char **args) 7 { 8 Mat A, B, C, At, Bt; 9 PetscViewer fd; 10 char file[PETSC_MAX_PATH_LEN]; 11 PetscRandom rand; 12 Vec xx, yy, s1, s2; 13 PetscReal s1norm, s2norm, rnorm, tol = 1.e-10; 14 PetscInt rstart, rend, rows[2], cols[2], m, n, i, j, M, N, ct, row, ncols1, ncols2, bs; 15 PetscMPIInt rank, size; 16 const PetscInt *cols1, *cols2; 17 PetscScalar vals1[4], vals2[4], v; 18 const PetscScalar *v1, *v2; 19 PetscBool flg; 20 21 PetscFunctionBeginUser; 22 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 23 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 24 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 25 26 /* Check out if MatLoad() works */ 27 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg)); 28 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Input file not specified"); 29 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 30 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 31 PetscCall(MatSetType(A, MATBAIJ)); 32 PetscCall(MatLoad(A, fd)); 33 34 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B)); 35 PetscCall(PetscViewerDestroy(&fd)); 36 37 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 38 PetscCall(PetscRandomSetFromOptions(rand)); 39 PetscCall(MatGetLocalSize(A, &m, &n)); 40 PetscCall(VecCreate(PETSC_COMM_WORLD, &xx)); 41 PetscCall(VecSetSizes(xx, m, PETSC_DECIDE)); 42 PetscCall(VecSetFromOptions(xx)); 43 PetscCall(VecDuplicate(xx, &s1)); 44 PetscCall(VecDuplicate(xx, &s2)); 45 PetscCall(VecDuplicate(xx, &yy)); 46 PetscCall(MatGetBlockSize(A, &bs)); 47 48 /* Test MatNorm() */ 49 PetscCall(MatNorm(A, NORM_FROBENIUS, &s1norm)); 50 PetscCall(MatNorm(B, NORM_FROBENIUS, &s2norm)); 51 rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm; 52 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs)); 53 PetscCall(MatNorm(A, NORM_INFINITY, &s1norm)); 54 PetscCall(MatNorm(B, NORM_INFINITY, &s2norm)); 55 rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm; 56 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs)); 57 PetscCall(MatNorm(A, NORM_1, &s1norm)); 58 PetscCall(MatNorm(B, NORM_1, &s2norm)); 59 rnorm = PetscAbsScalar(s2norm - s1norm) / s2norm; 60 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs)); 61 62 /* Test MatMult() */ 63 for (i = 0; i < IMAX; i++) { 64 PetscCall(VecSetRandom(xx, rand)); 65 PetscCall(MatMult(A, xx, s1)); 66 PetscCall(MatMult(B, xx, s2)); 67 PetscCall(VecAXPY(s2, -1.0, s1)); 68 PetscCall(VecNorm(s2, NORM_2, &rnorm)); 69 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMult - Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)rnorm, bs)); 70 } 71 72 /* test MatMultAdd() */ 73 for (i = 0; i < IMAX; i++) { 74 PetscCall(VecSetRandom(xx, rand)); 75 PetscCall(VecSetRandom(yy, rand)); 76 PetscCall(MatMultAdd(A, xx, yy, s1)); 77 PetscCall(MatMultAdd(B, xx, yy, s2)); 78 PetscCall(VecAXPY(s2, -1.0, s1)); 79 PetscCall(VecNorm(s2, NORM_2, &rnorm)); 80 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultAdd - Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)rnorm, bs)); 81 } 82 83 /* Test MatMultTranspose() */ 84 for (i = 0; i < IMAX; i++) { 85 PetscCall(VecSetRandom(xx, rand)); 86 PetscCall(MatMultTranspose(A, xx, s1)); 87 PetscCall(MatMultTranspose(B, xx, s2)); 88 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 89 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 90 rnorm = s2norm - s1norm; 91 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs)); 92 } 93 /* Test MatMultTransposeAdd() */ 94 for (i = 0; i < IMAX; i++) { 95 PetscCall(VecSetRandom(xx, rand)); 96 PetscCall(VecSetRandom(yy, rand)); 97 PetscCall(MatMultTransposeAdd(A, xx, yy, s1)); 98 PetscCall(MatMultTransposeAdd(B, xx, yy, s2)); 99 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 100 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 101 rnorm = s2norm - s1norm; 102 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs)); 103 } 104 105 /* Check MatGetValues() */ 106 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 107 PetscCall(MatGetSize(A, &M, &N)); 108 109 for (i = 0; i < IMAX; i++) { 110 /* Create random row numbers ad col numbers */ 111 PetscCall(PetscRandomGetValue(rand, &v)); 112 cols[0] = (int)(PetscRealPart(v) * N); 113 PetscCall(PetscRandomGetValue(rand, &v)); 114 cols[1] = (int)(PetscRealPart(v) * N); 115 PetscCall(PetscRandomGetValue(rand, &v)); 116 rows[0] = rstart + (int)(PetscRealPart(v) * m); 117 PetscCall(PetscRandomGetValue(rand, &v)); 118 rows[1] = rstart + (int)(PetscRealPart(v) * m); 119 120 PetscCall(MatGetValues(A, 2, rows, 2, cols, vals1)); 121 PetscCall(MatGetValues(B, 2, rows, 2, cols, vals2)); 122 123 for (j = 0; j < 4; j++) { 124 if (vals1[j] != vals2[j]) { 125 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]: Error: MatGetValues rstart = %2" PetscInt_FMT " row = %2" PetscInt_FMT " col = %2" PetscInt_FMT " val1 = %e val2 = %e bs = %" PetscInt_FMT "\n", rank, rstart, rows[j / 2], cols[j % 2], (double)PetscRealPart(vals1[j]), (double)PetscRealPart(vals2[j]), bs)); 126 } 127 } 128 } 129 130 /* Test MatGetRow()/ MatRestoreRow() */ 131 for (ct = 0; ct < 100; ct++) { 132 PetscCall(PetscRandomGetValue(rand, &v)); 133 row = rstart + (PetscInt)(PetscRealPart(v) * m); 134 PetscCall(MatGetRow(A, row, &ncols1, &cols1, &v1)); 135 PetscCall(MatGetRow(B, row, &ncols2, &cols2, &v2)); 136 137 for (i = 0, j = 0; i < ncols1 && j < ncols2; j++) { 138 while (cols2[j] != cols1[i]) i++; 139 PetscCheck(v1[i] == v2[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - vals incorrect."); 140 } 141 PetscCheck(j >= ncols2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - cols incorrect"); 142 143 PetscCall(MatRestoreRow(A, row, &ncols1, &cols1, &v1)); 144 PetscCall(MatRestoreRow(B, row, &ncols2, &cols2, &v2)); 145 } 146 147 /* Test MatConvert() */ 148 PetscCall(MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C)); 149 150 /* See if MatMult Says both are same */ 151 for (i = 0; i < IMAX; i++) { 152 PetscCall(VecSetRandom(xx, rand)); 153 PetscCall(MatMult(A, xx, s1)); 154 PetscCall(MatMult(C, xx, s2)); 155 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 156 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 157 rnorm = s2norm - s1norm; 158 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs)); 159 } 160 PetscCall(MatDestroy(&C)); 161 162 /* Test MatTranspose() */ 163 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At)); 164 PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &Bt)); 165 for (i = 0; i < IMAX; i++) { 166 PetscCall(VecSetRandom(xx, rand)); 167 PetscCall(MatMult(At, xx, s1)); 168 PetscCall(MatMult(Bt, xx, s2)); 169 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 170 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 171 rnorm = s2norm - s1norm; 172 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", rank, (double)s1norm, (double)s2norm, bs)); 173 } 174 PetscCall(MatDestroy(&At)); 175 PetscCall(MatDestroy(&Bt)); 176 177 PetscCall(MatDestroy(&A)); 178 PetscCall(MatDestroy(&B)); 179 PetscCall(VecDestroy(&xx)); 180 PetscCall(VecDestroy(&yy)); 181 PetscCall(VecDestroy(&s1)); 182 PetscCall(VecDestroy(&s2)); 183 PetscCall(PetscRandomDestroy(&rand)); 184 PetscCall(PetscFinalize()); 185 return 0; 186 } 187 188 /*TEST 189 190 build: 191 requires: !complex 192 193 test: 194 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 195 nsize: 3 196 args: -matload_block_size 1 -f ${DATAFILESPATH}/matrices/small 197 198 test: 199 suffix: 2 200 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 201 nsize: 3 202 args: -matload_block_size 2 -f ${DATAFILESPATH}/matrices/small 203 output_file: output/ex53_1.out 204 205 test: 206 suffix: 3 207 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 208 nsize: 3 209 args: -matload_block_size 4 -f ${DATAFILESPATH}/matrices/small 210 output_file: output/ex53_1.out 211 212 test: 213 TODO: Matrix row/column sizes are not compatible with block size 214 suffix: 4 215 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 216 nsize: 3 217 args: -matload_block_size 5 -f ${DATAFILESPATH}/matrices/small 218 output_file: output/ex53_1.out 219 220 test: 221 suffix: 5 222 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 223 nsize: 3 224 args: -matload_block_size 6 -f ${DATAFILESPATH}/matrices/small 225 output_file: output/ex53_1.out 226 227 test: 228 TODO: Matrix row/column sizes are not compatible with block size 229 suffix: 6 230 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 231 nsize: 3 232 args: -matload_block_size 7 -f ${DATAFILESPATH}/matrices/small 233 output_file: output/ex53_1.out 234 235 test: 236 TODO: Matrix row/column sizes are not compatible with block size 237 suffix: 7 238 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 239 nsize: 3 240 args: -matload_block_size 8 -f ${DATAFILESPATH}/matrices/small 241 output_file: output/ex53_1.out 242 243 test: 244 suffix: 8 245 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 246 nsize: 4 247 args: -matload_block_size 3 -f ${DATAFILESPATH}/matrices/small 248 output_file: output/ex53_1.out 249 250 TEST*/ 251