1 static char help[] = "Tests the use of MatTranspose_Nest and MatMatMult_Nest_Dense\n"; 2 3 #include <petscmat.h> 4 5 PetscErrorCode TestInitialMatrix(void) 6 { 7 const PetscInt nr = 2, nc = 3, nk = 10; 8 PetscInt n, N, m, M; 9 const PetscInt arow[2 * 3] = {2, 2, 2, 3, 3, 3}; 10 const PetscInt acol[2 * 3] = {3, 2, 4, 3, 2, 4}; 11 Mat A, Atranspose, B, C; 12 Mat subs[2 * 3], **block; 13 Vec x, y, Ax, ATy; 14 PetscInt i, j; 15 PetscScalar dot1, dot2, zero = 0.0, one = 1.0, *valsB, *valsC; 16 PetscReal norm; 17 PetscRandom rctx; 18 PetscBool equal; 19 20 PetscFunctionBegin; 21 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 22 /* Force the random numbers to have imaginary part 0 so printed results are the same for --with-scalar-type=real or --with-scalar-type=complex */ 23 PetscCall(PetscRandomSetInterval(rctx, zero, one)); 24 PetscCall(PetscRandomSetFromOptions(rctx)); 25 for (i = 0; i < (nr * nc); i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, arow[i], acol[i], NULL, &subs[i])); 26 PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, NULL, nc, NULL, subs, &A)); 27 PetscCall(MatCreateVecs(A, &x, NULL)); 28 PetscCall(MatCreateVecs(A, NULL, &y)); 29 PetscCall(VecDuplicate(x, &ATy)); 30 PetscCall(VecDuplicate(y, &Ax)); 31 PetscCall(MatSetRandom(A, rctx)); 32 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose)); 33 34 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 35 PetscCall(MatNestGetSubMats(A, NULL, NULL, &block)); 36 for (i = 0; i < nr; i++) { 37 for (j = 0; j < nc; j++) PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 38 } 39 40 PetscCall(MatView(Atranspose, PETSC_VIEWER_STDOUT_WORLD)); 41 PetscCall(MatNestGetSubMats(Atranspose, NULL, NULL, &block)); 42 for (i = 0; i < nc; i++) { 43 for (j = 0; j < nr; j++) PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 44 } 45 46 /* Check <Ax, y> = <x, A^Ty> */ 47 for (i = 0; i < 10; i++) { 48 PetscCall(VecSetRandom(x, rctx)); 49 PetscCall(VecSetRandom(y, rctx)); 50 51 PetscCall(MatMult(A, x, Ax)); 52 PetscCall(VecDot(Ax, y, &dot1)); 53 PetscCall(MatMult(Atranspose, y, ATy)); 54 PetscCall(VecDot(ATy, x, &dot2)); 55 56 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "<Ax, y> = %g\n", (double)PetscRealPart(dot1))); 57 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "<x, A^Ty> = %g\n", (double)PetscRealPart(dot2))); 58 } 59 PetscCall(VecDestroy(&x)); 60 PetscCall(VecDestroy(&y)); 61 PetscCall(VecDestroy(&Ax)); 62 63 PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, acol[0] + acol[nr] + acol[2 * nr], nk, NULL, &B)); 64 PetscCall(MatSetRandom(B, rctx)); 65 PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 66 PetscCall(MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C)); 67 PetscCall(MatMatMultEqual(A, B, C, 10, &equal)); 68 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in C != A*B"); 69 70 PetscCall(MatGetSize(A, &M, &N)); 71 PetscCall(MatGetLocalSize(A, &m, &n)); 72 for (i = 0; i < nk; i++) { 73 PetscCall(MatDenseGetColumn(B, i, &valsB)); 74 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, valsB, &x)); 75 PetscCall(MatCreateVecs(A, NULL, &Ax)); 76 PetscCall(MatMult(A, x, Ax)); 77 PetscCall(VecDestroy(&x)); 78 PetscCall(MatDenseGetColumn(C, i, &valsC)); 79 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, m, M, valsC, &y)); 80 PetscCall(VecAXPY(y, -1.0, Ax)); 81 PetscCall(VecDestroy(&Ax)); 82 PetscCall(VecDestroy(&y)); 83 PetscCall(MatDenseRestoreColumn(C, &valsC)); 84 PetscCall(MatDenseRestoreColumn(B, &valsB)); 85 } 86 PetscCall(MatNorm(C, NORM_INFINITY, &norm)); 87 PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatMatMult(): %g", (double)norm); 88 PetscCall(MatDestroy(&C)); 89 PetscCall(MatDestroy(&B)); 90 91 for (i = 0; i < (nr * nc); i++) PetscCall(MatDestroy(&subs[i])); 92 PetscCall(MatDestroy(&A)); 93 PetscCall(MatDestroy(&Atranspose)); 94 PetscCall(VecDestroy(&ATy)); 95 PetscCall(PetscRandomDestroy(&rctx)); 96 PetscFunctionReturn(0); 97 } 98 99 PetscErrorCode TestReuseMatrix(void) 100 { 101 const PetscInt n = 2; 102 Mat A; 103 Mat subs[2 * 2], **block; 104 PetscInt i, j; 105 PetscRandom rctx; 106 PetscScalar zero = 0.0, one = 1.0; 107 108 PetscFunctionBegin; 109 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 110 PetscCall(PetscRandomSetInterval(rctx, zero, one)); 111 PetscCall(PetscRandomSetFromOptions(rctx)); 112 for (i = 0; i < (n * n); i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, n, n, NULL, &subs[i])); 113 PetscCall(MatCreateNest(PETSC_COMM_WORLD, n, NULL, n, NULL, subs, &A)); 114 PetscCall(MatSetRandom(A, rctx)); 115 116 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 117 PetscCall(MatNestGetSubMats(A, NULL, NULL, &block)); 118 for (i = 0; i < n; i++) { 119 for (j = 0; j < n; j++) PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 120 } 121 PetscCall(MatTranspose(A, MAT_INPLACE_MATRIX, &A)); 122 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 123 PetscCall(MatNestGetSubMats(A, NULL, NULL, &block)); 124 for (i = 0; i < n; i++) { 125 for (j = 0; j < n; j++) PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 126 } 127 128 for (i = 0; i < (n * n); i++) PetscCall(MatDestroy(&subs[i])); 129 PetscCall(MatDestroy(&A)); 130 PetscCall(PetscRandomDestroy(&rctx)); 131 PetscFunctionReturn(0); 132 } 133 134 int main(int argc, char **args) 135 { 136 PetscFunctionBeginUser; 137 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 138 PetscCall(TestInitialMatrix()); 139 PetscCall(TestReuseMatrix()); 140 PetscCall(PetscFinalize()); 141 return 0; 142 } 143 144 /*TEST 145 146 test: 147 args: -malloc_dump 148 149 TEST*/ 150