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