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++) {
25: MatCreateSeqDense(PETSC_COMM_WORLD,arow[i],acol[i],NULL,&subs[i]);
26: }
27: MatCreateNest(PETSC_COMM_WORLD,nr,NULL,nc,NULL,subs,&A);
28: MatCreateVecs(A, &x, NULL);
29: MatCreateVecs(A, NULL, &y);
30: VecDuplicate(x, &ATy);
31: VecDuplicate(y, &Ax);
32: MatSetRandom(A,rctx);
33: MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose);
35: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
36: MatNestGetSubMats(A, NULL, NULL, &block);
37: for (i=0; i<nr; i++) {
38: for (j=0; j<nc; j++) {
39: MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
40: }
41: }
43: MatView(Atranspose, PETSC_VIEWER_STDOUT_WORLD);
44: MatNestGetSubMats(Atranspose, NULL, NULL, &block);
45: for (i=0; i<nc; i++) {
46: for (j=0; j<nr; j++) {
47: MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
48: }
49: }
51: /* Check <Ax, y> = <x, A^Ty> */
52: for (i=0; i<10; i++) {
53: VecSetRandom(x,rctx);
54: VecSetRandom(y,rctx);
56: MatMult(A, x, Ax);
57: VecDot(Ax, y, &dot1);
58: MatMult(Atranspose, y, ATy);
59: VecDot(ATy, x, &dot2);
61: PetscPrintf(PETSC_COMM_WORLD, "<Ax, y> = %g\n", (double)PetscRealPart(dot1));
62: PetscPrintf(PETSC_COMM_WORLD, "<x, A^Ty> = %g\n",(double)PetscRealPart(dot2));
63: }
64: VecDestroy(&x);
65: VecDestroy(&y);
66: VecDestroy(&Ax);
68: MatCreateSeqDense(PETSC_COMM_WORLD,acol[0]+acol[nr]+acol[2*nr],nk,NULL,&B);
69: MatSetRandom(B,rctx);
70: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
71: MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
72: MatMatMultEqual(A,B,C,10,&equal);
75: MatGetSize(A,&M,&N);
76: MatGetLocalSize(A,&m,&n);
77: for (i=0; i<nk; i++) {
78: MatDenseGetColumn(B,i,&valsB);
79: VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,valsB,&x);
80: MatCreateVecs(A,NULL,&Ax);
81: MatMult(A,x,Ax);
82: VecDestroy(&x);
83: MatDenseGetColumn(C,i,&valsC);
84: VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m,M,valsC,&y);
85: VecAXPY(y,-1.0,Ax);
86: VecDestroy(&Ax);
87: VecDestroy(&y);
88: MatDenseRestoreColumn(C,&valsC);
89: MatDenseRestoreColumn(B,&valsB);
90: }
91: MatNorm(C,NORM_INFINITY,&norm);
93: MatDestroy(&C);
94: MatDestroy(&B);
96: for (i=0; i<(nr * nc); i++) {
97: MatDestroy(&subs[i]);
98: }
99: MatDestroy(&A);
100: MatDestroy(&Atranspose);
101: VecDestroy(&ATy);
102: PetscRandomDestroy(&rctx);
103: return 0;
104: }
106: PetscErrorCode TestReuseMatrix(void)
107: {
108: const PetscInt n = 2;
109: Mat A;
110: Mat subs[2*2],**block;
111: PetscInt i,j;
112: PetscRandom rctx;
113: PetscScalar zero = 0.0, one = 1.0;
115: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
116: PetscRandomSetInterval(rctx,zero,one);
117: PetscRandomSetFromOptions(rctx);
118: for (i=0; i<(n * n); i++) {
119: MatCreateSeqDense(PETSC_COMM_WORLD,n,n,NULL,&subs[i]);
120: }
121: MatCreateNest(PETSC_COMM_WORLD,n,NULL,n,NULL,subs,&A);
122: MatSetRandom(A,rctx);
124: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
125: MatNestGetSubMats(A, NULL, NULL, &block);
126: for (i=0; i<n; i++) {
127: for (j=0; j<n; j++) {
128: MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
129: }
130: }
131: MatTranspose(A,MAT_INPLACE_MATRIX,&A);
132: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
133: MatNestGetSubMats(A, NULL, NULL, &block);
134: for (i=0; i<n; i++) {
135: for (j=0; j<n; j++) {
136: MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD);
137: }
138: }
140: for (i=0; i<(n * n); i++) {
141: MatDestroy(&subs[i]);
142: }
143: MatDestroy(&A);
144: PetscRandomDestroy(&rctx);
145: return 0;
146: }
148: int main(int argc,char **args)
149: {
151: PetscInitialize(&argc,&args,(char*)0,help);
152: TestInitialMatrix();
153: TestReuseMatrix();
154: PetscFinalize();
155: return 0;
156: }
158: /*TEST
160: test:
161: args: -malloc_dump
163: TEST*/