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