1 static char help[] = "Test MatMatMult(), MatTranspose(), MatTransposeMatMult() for Dense and Elemental matrices.\n\n"; 2 /* 3 Example: 4 mpiexec -n <np> ./ex104 -mat_type elemental 5 */ 6 7 #include <petscmat.h> 8 9 int main(int argc, char **argv) 10 { 11 Mat A, B, C, D; 12 PetscInt i, M = 10, N = 5, j, nrows, ncols, am, an, rstart, rend; 13 PetscRandom r; 14 PetscBool equal, Aiselemental; 15 PetscBool columns_on_one_rank = PETSC_FALSE; 16 PetscReal fill = 1.0; 17 IS isrows, iscols; 18 const PetscInt *rows, *cols; 19 PetscScalar *v, rval; 20 PetscBool Test_MatMatMult = PETSC_TRUE; 21 PetscMPIInt size, rank; 22 23 PetscFunctionBeginUser; 24 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 25 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 26 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 27 28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 30 PetscCall(PetscOptionsGetBool(NULL, NULL, "-columns_on_one_rank", &columns_on_one_rank, NULL)); 31 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 32 if (!columns_on_one_rank) { 33 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N)); 34 } else { 35 PetscCall(MatSetSizes(A, PETSC_DECIDE, rank == 0 ? N : 0, M, N)); 36 } 37 PetscCall(MatSetType(A, MATDENSE)); 38 PetscCall(MatSetFromOptions(A)); 39 PetscCall(MatSetUp(A)); 40 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r)); 41 PetscCall(PetscRandomSetFromOptions(r)); 42 43 /* Set local matrix entries */ 44 PetscCall(MatGetOwnershipIS(A, &isrows, &iscols)); 45 PetscCall(ISGetLocalSize(isrows, &nrows)); 46 PetscCall(ISGetIndices(isrows, &rows)); 47 PetscCall(ISGetLocalSize(iscols, &ncols)); 48 PetscCall(ISGetIndices(iscols, &cols)); 49 PetscCall(PetscMalloc1(nrows * ncols, &v)); 50 for (i = 0; i < nrows; i++) { 51 for (j = 0; j < ncols; j++) { 52 PetscCall(PetscRandomGetValue(r, &rval)); 53 v[i * ncols + j] = rval; 54 } 55 } 56 PetscCall(MatSetValues(A, nrows, rows, ncols, cols, v, INSERT_VALUES)); 57 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 58 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 59 PetscCall(ISRestoreIndices(isrows, &rows)); 60 PetscCall(ISRestoreIndices(iscols, &cols)); 61 PetscCall(ISDestroy(&isrows)); 62 PetscCall(ISDestroy(&iscols)); 63 PetscCall(PetscRandomDestroy(&r)); 64 65 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATELEMENTAL, &Aiselemental)); 66 67 /* Test MatCreateTranspose() and MatTranspose() */ 68 PetscCall(MatCreateTranspose(A, &C)); 69 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B)); /* B = A^T */ 70 PetscCall(MatMultEqual(C, B, 10, &equal)); 71 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A^T*x != (x^T*A)^T"); 72 PetscCall(MatDestroy(&B)); 73 74 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B)); 75 if (!Aiselemental) { 76 PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B)); 77 PetscCall(MatMultEqual(C, B, 10, &equal)); 78 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "C*x != B*x"); 79 } 80 PetscCall(MatDestroy(&B)); 81 82 /* Test B = C*A for matrix type transpose and seqdense */ 83 if (size == 1 && !Aiselemental) { 84 PetscCall(MatScale(C, -1.0)); 85 PetscCall(MatMatMult(C, A, MAT_INITIAL_MATRIX, fill, &B)); 86 PetscCall(MatScale(C, -1.0)); 87 PetscCall(MatScale(B, -1.0)); 88 PetscCall(MatMatMultEqual(C, A, B, 10, &equal)); 89 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "B != C*A for matrix type transpose and seqdense"); 90 PetscCall(MatDestroy(&B)); 91 } 92 PetscCall(MatDestroy(&C)); 93 94 /* Test MatMatMult() */ 95 if (Test_MatMatMult) { 96 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B)); /* B = A^T */ 97 PetscCall(MatMatMult(B, A, MAT_INITIAL_MATRIX, fill, &C)); /* C = B*A = A^T*A */ 98 PetscCall(MatMatMult(B, A, MAT_REUSE_MATRIX, fill, &C)); 99 PetscCall(MatMatMultEqual(B, A, C, 10, &equal)); 100 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "B*A*x != C*x"); 101 102 /* Test MatDuplicate for matrix product */ 103 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &D)); 104 105 PetscCall(MatDestroy(&D)); 106 PetscCall(MatDestroy(&C)); 107 PetscCall(MatDestroy(&B)); 108 } 109 110 /* Test MatTransposeMatMult() */ 111 if (!Aiselemental) { 112 Mat E; 113 114 PetscCall(MatTransposeMatMult(A, A, MAT_INITIAL_MATRIX, fill, &D)); /* D = A^T*A */ 115 PetscCall(MatTransposeMatMult(A, A, MAT_REUSE_MATRIX, fill, &D)); 116 PetscCall(MatTransposeMatMultEqual(A, A, D, 10, &equal)); 117 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "D*x != A^T*A*x"); 118 119 /* Test MatDuplicate for matrix product */ 120 PetscCall(MatDuplicate(D, MAT_COPY_VALUES, &C)); 121 PetscCall(MatDestroy(&C)); 122 123 /* Test A*D for fast path when D is on one process */ 124 PetscCall(MatSetRandom(D, NULL)); 125 PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, fill, &E)); 126 PetscCall(MatMatMult(A, D, MAT_REUSE_MATRIX, fill, &E)); 127 PetscCall(MatMatMultEqual(A, D, E, 10, &equal)); 128 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "E*x != A*D*x"); 129 PetscCall(MatDestroy(&E)); 130 131 PetscCall(MatDestroy(&D)); 132 133 /* Test D*x = A^T*C*A*x, where C is in AIJ format */ 134 PetscCall(MatGetLocalSize(A, &am, &an)); 135 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 136 if (size == 1) { 137 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, am, am)); 138 } else { 139 PetscCall(MatSetSizes(C, am, am, PETSC_DECIDE, PETSC_DECIDE)); 140 } 141 PetscCall(MatSetFromOptions(C)); 142 PetscCall(MatSetUp(C)); 143 PetscCall(MatGetOwnershipRange(C, &rstart, &rend)); 144 v[0] = 1.0; 145 for (i = rstart; i < rend; i++) PetscCall(MatSetValues(C, 1, &i, 1, &i, v, INSERT_VALUES)); 146 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 147 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 148 149 /* B = C*A, D = A^T*B */ 150 PetscCall(MatMatMult(C, A, MAT_INITIAL_MATRIX, 1.0, &B)); 151 PetscCall(MatTransposeMatMult(A, B, MAT_INITIAL_MATRIX, fill, &D)); 152 PetscCall(MatTransposeMatMultEqual(A, B, D, 10, &equal)); 153 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "D*x != A^T*B*x"); 154 155 PetscCall(MatDestroy(&D)); 156 PetscCall(MatDestroy(&C)); 157 PetscCall(MatDestroy(&B)); 158 } 159 160 /* Test MatMatTransposeMult() */ 161 if (!Aiselemental) { 162 PetscReal diff, scale; 163 PetscInt am, an, aM, aN; 164 165 PetscCall(MatGetLocalSize(A, &am, &an)); 166 PetscCall(MatGetSize(A, &aM, &aN)); 167 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), PETSC_DECIDE, an, aM + 10, aN, NULL, &B)); 168 PetscCall(MatSetRandom(B, NULL)); 169 PetscCall(MatMatTransposeMult(A, B, MAT_INITIAL_MATRIX, fill, &D)); /* D = A*A^T */ 170 171 /* Test MatDuplicate for matrix product */ 172 PetscCall(MatDuplicate(D, MAT_COPY_VALUES, &C)); 173 174 PetscCall(MatMatTransposeMult(A, B, MAT_REUSE_MATRIX, fill, &D)); 175 PetscCall(MatAXPY(C, -1., D, SAME_NONZERO_PATTERN)); 176 177 PetscCall(MatNorm(C, NORM_FROBENIUS, &diff)); 178 PetscCall(MatNorm(D, NORM_FROBENIUS, &scale)); 179 PetscCheck(diff <= PETSC_SMALL * scale, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatMatTransposeMult() differs between MAT_INITIAL_MATRIX and MAT_REUSE_MATRIX"); 180 PetscCall(MatDestroy(&C)); 181 182 PetscCall(MatMatTransposeMultEqual(A, B, D, 10, &equal)); 183 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "D*x != A^T*A*x"); 184 PetscCall(MatDestroy(&D)); 185 PetscCall(MatDestroy(&B)); 186 } 187 188 PetscCall(MatDestroy(&A)); 189 PetscCall(PetscFree(v)); 190 PetscCall(PetscFinalize()); 191 return 0; 192 } 193 194 /*TEST 195 196 test: 197 output_file: output/empty.out 198 199 test: 200 suffix: 2 201 nsize: 2 202 output_file: output/empty.out 203 204 test: 205 suffix: 3 206 nsize: 4 207 output_file: output/empty.out 208 args: -M 23 -N 31 209 210 test: 211 suffix: 4 212 nsize: 4 213 output_file: output/empty.out 214 args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via cyclic 215 216 test: 217 suffix: 5 218 nsize: 4 219 output_file: output/empty.out 220 args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via allgatherv 221 222 test: 223 suffix: 6 224 args: -mat_type elemental 225 requires: elemental 226 output_file: output/empty.out 227 228 testset: 229 nsize: 2 230 output_file: output/empty.out 231 requires: elemental 232 test: 233 suffix: 7_dense 234 args: -mat_type dense -mat_product_algorithm elemental 235 test: 236 suffix: 7_elemental 237 args: -mat_type elemental 238 239 test: 240 suffix: 8 241 nsize: 4 242 args: -columns_on_one_rank 243 output_file: output/empty.out 244 245 test: 246 suffix: 9 247 nsize: 4 248 requires: cuda 249 args: -columns_on_one_rank -mat_type densecuda 250 output_file: output/empty.out 251 252 TEST*/ 253