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;

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

 23:   /* Determine file from which we read the matrix A */
 24:   PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg);

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

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

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

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

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

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

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

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

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

111:     MatMultTranspose(A, x, y);

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

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

130:   /* Free data structures */
131:   MatDestroy(&A);
132:   MatDestroy(&C);
133:   MatDestroy(&Bdense);
134:   MatDestroy(&Cdense);
135:   VecDestroy(&x);
136:   VecDestroy(&y);
137:   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*/