Actual source code: ex202.c

  1: static char help[] = "Tests the use of MatTranspose_Nest and MatMatMult_Nest_Dense\n";

  3: #include <petscmat.h>

  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;

 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));

 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:   }

 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:   }

 46:   /* Check <Ax, y> = <x, A^Ty> */
 47:   for (i = 0; i < 10; i++) {
 48:     PetscCall(VecSetRandom(x, rctx));
 49:     PetscCall(VecSetRandom(y, rctx));

 51:     PetscCall(MatMult(A, x, Ax));
 52:     PetscCall(VecDot(Ax, y, &dot1));
 53:     PetscCall(MatMult(Atranspose, y, ATy));
 54:     PetscCall(VecDot(ATy, x, &dot2));

 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));

 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");

 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));

 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(PETSC_SUCCESS);
 97: }

 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;

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));

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:   }

128:   for (i = 0; i < (n * n); i++) PetscCall(MatDestroy(&subs[i]));
129:   PetscCall(MatDestroy(&A));
130:   PetscCall(PetscRandomDestroy(&rctx));
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

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: }

144: /*TEST

146:    test:
147:       args: -malloc_dump

149: TEST*/