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