Actual source code: ex209.c
1: static char help[] = "Test MatTransposeMatMult() \n\n";
3: /* Example:
4: mpiexec -n 8 ./ex209 -f <inputfile> (e.g., inputfile=ceres_solver_iteration_001_A.petsc)
5: */
7: #include <petscmat.h>
9: int main(int argc,char **args)
10: {
11: Mat A,C,AtA,B;
12: PetscViewer fd;
13: char file[PETSC_MAX_PATH_LEN];
14: PetscReal fill = 4.0;
15: PetscMPIInt rank,size;
16: PetscBool equal;
17: PetscInt i,m,n,rstart,rend;
18: PetscScalar one=1.0;
20: PetscInitialize(&argc,&args,(char*)0,help);
21: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
22: MPI_Comm_size(PETSC_COMM_WORLD,&size);
23: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),NULL);
25: /* Load matrix A */
26: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
27: MatCreate(PETSC_COMM_WORLD,&A);
28: MatSetType(A,MATAIJ);
29: MatLoad(A,fd);
30: PetscViewerDestroy(&fd);
32: /* Create identity matrix B */
33: MatGetLocalSize(A,&m,&n);
34: MatCreate(PETSC_COMM_WORLD,&B);
35: MatSetSizes(B,m,m,PETSC_DECIDE,PETSC_DECIDE);
36: MatSetType(B,MATAIJ);
37: MatSetFromOptions(B);
38: MatSetUp(B);
40: MatGetOwnershipRange(B,&rstart,&rend);
41: for (i=rstart; i<rend; i++) {
42: MatSetValues(B,1,&i,1,&i,&one,INSERT_VALUES);
43: }
44: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
45: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
47: /* Compute AtA = A^T*B*A, B = identity matrix */
48: MatPtAP(B,A,MAT_INITIAL_MATRIX,fill,&AtA);
49: MatPtAP(B,A,MAT_REUSE_MATRIX,fill,&AtA);
50: if (rank == 0) printf("C = A^T*B*A is done...\n");
51: MatDestroy(&B);
53: /* Compute C = A^T*A */
54: MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&C);
55: if (rank == 0) printf("C = A^T*A is done...\n");
56: MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&C);
57: if (rank == 0) printf("REUSE C = A^T*A is done...\n");
59: /* Compare C and AtA */
60: MatMultEqual(C,AtA,20,&equal);
62: MatDestroy(&AtA);
64: MatDestroy(&C);
65: MatDestroy(&A);
66: PetscFinalize();
67: return 0;
68: }
70: /*TEST
72: test:
73: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
74: args: -f ${DATAFILESPATH}/matrices/arco1
76: test:
77: suffix: 2
78: nsize: 4
79: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
80: args: -f ${DATAFILESPATH}/matrices/arco1 -matptap_via nonscalable
82: TEST*/