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:   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:   PetscRandomSetInterval(rctx, zero, one);
 23:   PetscRandomSetFromOptions(rctx);
 24:   for (i = 0; i < (nr * nc); i++) MatCreateSeqDense(PETSC_COMM_WORLD, arow[i], acol[i], NULL, &subs[i]);
 25:   MatCreateNest(PETSC_COMM_WORLD, nr, NULL, nc, NULL, subs, &A);
 26:   MatCreateVecs(A, &x, NULL);
 27:   MatCreateVecs(A, NULL, &y);
 28:   VecDuplicate(x, &ATy);
 29:   VecDuplicate(y, &Ax);
 30:   MatSetRandom(A, rctx);
 31:   MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose);

 33:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
 34:   MatNestGetSubMats(A, NULL, NULL, &block);
 35:   for (i = 0; i < nr; i++) {
 36:     for (j = 0; j < nc; j++) MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
 37:   }

 39:   MatView(Atranspose, PETSC_VIEWER_STDOUT_WORLD);
 40:   MatNestGetSubMats(Atranspose, NULL, NULL, &block);
 41:   for (i = 0; i < nc; i++) {
 42:     for (j = 0; j < nr; j++) MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
 43:   }

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

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

 55:     PetscPrintf(PETSC_COMM_WORLD, "<Ax, y> = %g\n", (double)PetscRealPart(dot1));
 56:     PetscPrintf(PETSC_COMM_WORLD, "<x, A^Ty> = %g\n", (double)PetscRealPart(dot2));
 57:   }
 58:   VecDestroy(&x);
 59:   VecDestroy(&y);
 60:   VecDestroy(&Ax);

 62:   MatCreateSeqDense(PETSC_COMM_WORLD, acol[0] + acol[nr] + acol[2 * nr], nk, NULL, &B);
 63:   MatSetRandom(B, rctx);
 64:   MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
 65:   MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);
 66:   MatMatMultEqual(A, B, C, 10, &equal);

 69:   MatGetSize(A, &M, &N);
 70:   MatGetLocalSize(A, &m, &n);
 71:   for (i = 0; i < nk; i++) {
 72:     MatDenseGetColumn(B, i, &valsB);
 73:     VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, valsB, &x);
 74:     MatCreateVecs(A, NULL, &Ax);
 75:     MatMult(A, x, Ax);
 76:     VecDestroy(&x);
 77:     MatDenseGetColumn(C, i, &valsC);
 78:     VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, m, M, valsC, &y);
 79:     VecAXPY(y, -1.0, Ax);
 80:     VecDestroy(&Ax);
 81:     VecDestroy(&y);
 82:     MatDenseRestoreColumn(C, &valsC);
 83:     MatDenseRestoreColumn(B, &valsB);
 84:   }
 85:   MatNorm(C, NORM_INFINITY, &norm);
 87:   MatDestroy(&C);
 88:   MatDestroy(&B);

 90:   for (i = 0; i < (nr * nc); i++) MatDestroy(&subs[i]);
 91:   MatDestroy(&A);
 92:   MatDestroy(&Atranspose);
 93:   VecDestroy(&ATy);
 94:   PetscRandomDestroy(&rctx);
 95:   return 0;
 96: }

 98: PetscErrorCode TestReuseMatrix(void)
 99: {
100:   const PetscInt n = 2;
101:   Mat            A;
102:   Mat            subs[2 * 2], **block;
103:   PetscInt       i, j;
104:   PetscRandom    rctx;
105:   PetscScalar    zero = 0.0, one = 1.0;

107:   PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
108:   PetscRandomSetInterval(rctx, zero, one);
109:   PetscRandomSetFromOptions(rctx);
110:   for (i = 0; i < (n * n); i++) MatCreateSeqDense(PETSC_COMM_WORLD, n, n, NULL, &subs[i]);
111:   MatCreateNest(PETSC_COMM_WORLD, n, NULL, n, NULL, subs, &A);
112:   MatSetRandom(A, rctx);

114:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
115:   MatNestGetSubMats(A, NULL, NULL, &block);
116:   for (i = 0; i < n; i++) {
117:     for (j = 0; j < n; j++) MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
118:   }
119:   MatTranspose(A, MAT_INPLACE_MATRIX, &A);
120:   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
121:   MatNestGetSubMats(A, NULL, NULL, &block);
122:   for (i = 0; i < n; i++) {
123:     for (j = 0; j < n; j++) MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
124:   }

126:   for (i = 0; i < (n * n); i++) MatDestroy(&subs[i]);
127:   MatDestroy(&A);
128:   PetscRandomDestroy(&rctx);
129:   return 0;
130: }

132: int main(int argc, char **args)
133: {
135:   PetscInitialize(&argc, &args, (char *)0, help);
136:   TestInitialMatrix();
137:   TestReuseMatrix();
138:   PetscFinalize();
139:   return 0;
140: }

142: /*TEST

144:    test:
145:       args: -malloc_dump

147: TEST*/