Actual source code: ex163.c


  2: static char help[] = "Tests MatTransposeMatMult() on MatLoad() matrix \n\n";

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Mat          A, C, Bdense, Cdense;
  9:   PetscViewer  fd;                       /* viewer */
 10:   char         file[PETSC_MAX_PATH_LEN]; /* input file name */
 11:   PetscBool    flg, viewmats = PETSC_FALSE;
 12:   PetscMPIInt  rank, size;
 13:   PetscReal    fill = 1.0;
 14:   PetscInt     m, n, i, j, BN = 10, rstart, rend, *rows, *cols;
 15:   PetscScalar *Barray, *Carray, rval, *array;
 16:   Vec          x, y;
 17:   PetscRandom  rand;

 19:   PetscFunctionBeginUser;
 20:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 21:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 23:   /* Determine file from which we read the matrix A */
 24:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
 25:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");

 27:   /* Load matrix A */
 28:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 29:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 30:   PetscCall(MatLoad(A, fd));
 31:   PetscCall(PetscViewerDestroy(&fd));

 33:   /* Print (for testing only) */
 34:   PetscCall(PetscOptionsHasName(NULL, NULL, "-view_mats", &viewmats));
 35:   if (viewmats) {
 36:     if (rank == 0) printf("A_aij:\n");
 37:     PetscCall(MatView(A, 0));
 38:   }

 40:   /* Test MatTransposeMatMult_aij_aij() */
 41:   PetscCall(MatTransposeMatMult(A, A, MAT_INITIAL_MATRIX, fill, &C));
 42:   if (viewmats) {
 43:     if (rank == 0) printf("\nC = A_aij^T * A_aij:\n");
 44:     PetscCall(MatView(C, 0));
 45:   }
 46:   PetscCall(MatDestroy(&C));
 47:   PetscCall(MatGetLocalSize(A, &m, &n));

 49:   /* create a dense matrix Bdense */
 50:   PetscCall(MatCreate(PETSC_COMM_WORLD, &Bdense));
 51:   PetscCall(MatSetSizes(Bdense, m, PETSC_DECIDE, PETSC_DECIDE, BN));
 52:   PetscCall(MatSetType(Bdense, MATDENSE));
 53:   PetscCall(MatSetFromOptions(Bdense));
 54:   PetscCall(MatSetUp(Bdense));
 55:   PetscCall(MatGetOwnershipRange(Bdense, &rstart, &rend));

 57:   PetscCall(PetscMalloc3(m, &rows, BN, &cols, m * BN, &array));
 58:   for (i = 0; i < m; i++) rows[i] = rstart + i;
 59:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 60:   PetscCall(PetscRandomSetFromOptions(rand));
 61:   for (j = 0; j < BN; j++) {
 62:     cols[j] = j;
 63:     for (i = 0; i < m; i++) {
 64:       PetscCall(PetscRandomGetValue(rand, &rval));
 65:       array[m * j + i] = rval;
 66:     }
 67:   }
 68:   PetscCall(MatSetValues(Bdense, m, rows, BN, cols, array, INSERT_VALUES));
 69:   PetscCall(MatAssemblyBegin(Bdense, MAT_FINAL_ASSEMBLY));
 70:   PetscCall(MatAssemblyEnd(Bdense, MAT_FINAL_ASSEMBLY));
 71:   PetscCall(PetscRandomDestroy(&rand));
 72:   PetscCall(PetscFree3(rows, cols, array));
 73:   if (viewmats) {
 74:     if (rank == 0) printf("\nBdense:\n");
 75:     PetscCall(MatView(Bdense, 0));
 76:   }

 78:   /* Test MatTransposeMatMult_aij_dense() */
 79:   PetscCall(MatTransposeMatMult(A, Bdense, MAT_INITIAL_MATRIX, fill, &C));
 80:   PetscCall(MatTransposeMatMult(A, Bdense, MAT_REUSE_MATRIX, fill, &C));
 81:   if (viewmats) {
 82:     if (rank == 0) printf("\nC=A^T*Bdense:\n");
 83:     PetscCall(MatView(C, 0));
 84:   }

 86:   /* Check accuracy */
 87:   PetscCall(MatCreate(PETSC_COMM_WORLD, &Cdense));
 88:   PetscCall(MatSetSizes(Cdense, n, PETSC_DECIDE, PETSC_DECIDE, BN));
 89:   PetscCall(MatSetType(Cdense, MATDENSE));
 90:   PetscCall(MatSetFromOptions(Cdense));
 91:   PetscCall(MatSetUp(Cdense));
 92:   PetscCall(MatAssemblyBegin(Cdense, MAT_FINAL_ASSEMBLY));
 93:   PetscCall(MatAssemblyEnd(Cdense, MAT_FINAL_ASSEMBLY));

 95:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 96:   if (size == 1) {
 97:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &x));
 98:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n, NULL, &y));
 99:   } else {
100:     PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, m, PETSC_DECIDE, NULL, &x));
101:     PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, PETSC_DECIDE, NULL, &y));
102:   }

104:   /* Cdense[:,j] = A^T * Bdense[:,j] */
105:   PetscCall(MatDenseGetArray(Bdense, &Barray));
106:   PetscCall(MatDenseGetArray(Cdense, &Carray));
107:   for (j = 0; j < BN; j++) {
108:     PetscCall(VecPlaceArray(x, Barray));
109:     PetscCall(VecPlaceArray(y, Carray));

111:     PetscCall(MatMultTranspose(A, x, y));

113:     PetscCall(VecResetArray(x));
114:     PetscCall(VecResetArray(y));
115:     Barray += m;
116:     Carray += n;
117:   }
118:   PetscCall(MatDenseRestoreArray(Bdense, &Barray));
119:   PetscCall(MatDenseRestoreArray(Cdense, &Carray));
120:   if (viewmats) {
121:     if (rank == 0) printf("\nCdense:\n");
122:     PetscCall(MatView(Cdense, 0));
123:   }

125:   PetscCall(MatEqual(C, Cdense, &flg));
126:   if (!flg) {
127:     if (rank == 0) printf(" C != Cdense\n");
128:   }

130:   /* Free data structures */
131:   PetscCall(MatDestroy(&A));
132:   PetscCall(MatDestroy(&C));
133:   PetscCall(MatDestroy(&Bdense));
134:   PetscCall(MatDestroy(&Cdense));
135:   PetscCall(VecDestroy(&x));
136:   PetscCall(VecDestroy(&y));
137:   PetscCall(PetscFinalize());
138:   return 0;
139: }

141: /*TEST

143:    test:
144:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
145:       args: -f ${DATAFILESPATH}/matrices/small
146:       output_file: output/ex163.out

148:    test:
149:       suffix: 2
150:       nsize: 3
151:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
152:       args: -f ${DATAFILESPATH}/matrices/small
153:       output_file: output/ex163.out

155: TEST*/