Actual source code: ex104.c
1: static char help[] = "Test MatMatMult(), MatTranspose(), MatTransposeMatMult() for Dense and Elemental matrices.\n\n";
2: /*
3: Example:
4: mpiexec -n <np> ./ex104 -mat_type elemental
5: */
7: #include <petscmat.h>
9: int main(int argc,char **argv)
10: {
11: Mat A,B,C,D;
12: PetscInt i,M=10,N=5,j,nrows,ncols,am,an,rstart,rend;
13: PetscRandom r;
14: PetscBool equal,Aiselemental;
15: PetscReal fill = 1.0;
16: IS isrows,iscols;
17: const PetscInt *rows,*cols;
18: PetscScalar *v,rval;
19: #if defined(PETSC_HAVE_ELEMENTAL)
20: PetscBool Test_MatMatMult=PETSC_TRUE;
21: #else
22: PetscBool Test_MatMatMult=PETSC_FALSE;
23: #endif
24: PetscMPIInt size;
26: PetscInitialize(&argc,&argv,(char*)0,help);
27: MPI_Comm_size(PETSC_COMM_WORLD,&size);
29: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
30: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
31: MatCreate(PETSC_COMM_WORLD,&A);
32: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
33: MatSetType(A,MATDENSE);
34: MatSetFromOptions(A);
35: MatSetUp(A);
36: PetscRandomCreate(PETSC_COMM_WORLD,&r);
37: PetscRandomSetFromOptions(r);
39: /* Set local matrix entries */
40: MatGetOwnershipIS(A,&isrows,&iscols);
41: ISGetLocalSize(isrows,&nrows);
42: ISGetIndices(isrows,&rows);
43: ISGetLocalSize(iscols,&ncols);
44: ISGetIndices(iscols,&cols);
45: PetscMalloc1(nrows*ncols,&v);
46: for (i=0; i<nrows; i++) {
47: for (j=0; j<ncols; j++) {
48: PetscRandomGetValue(r,&rval);
49: v[i*ncols+j] = rval;
50: }
51: }
52: MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES);
53: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
54: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
55: ISRestoreIndices(isrows,&rows);
56: ISRestoreIndices(iscols,&cols);
57: ISDestroy(&isrows);
58: ISDestroy(&iscols);
59: PetscRandomDestroy(&r);
61: PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&Aiselemental);
63: /* Test MatCreateTranspose() and MatTranspose() */
64: MatCreateTranspose(A,&C);
65: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
66: MatMultEqual(C,B,10,&equal);
68: MatDestroy(&B);
70: MatDuplicate(A,MAT_COPY_VALUES,&B);
71: if (!Aiselemental) {
72: MatTranspose(B,MAT_INPLACE_MATRIX,&B);
73: MatMultEqual(C,B,10,&equal);
75: }
76: MatDestroy(&B);
78: /* Test B = C*A for matrix type transpose and seqdense */
79: if (size == 1 && !Aiselemental) {
80: MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&B);
81: MatMatMultEqual(C,A,B,10,&equal);
83: MatDestroy(&B);
84: }
85: MatDestroy(&C);
87: /* Test MatMatMult() */
88: if (Test_MatMatMult) {
89: #if !defined(PETSC_HAVE_ELEMENTAL)
91: #endif
92: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
93: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */
94: MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C);
95: MatMatMultEqual(B,A,C,10,&equal);
98: /* Test MatDuplicate for matrix product */
99: MatDuplicate(C,MAT_COPY_VALUES,&D);
101: MatDestroy(&D);
102: MatDestroy(&C);
103: MatDestroy(&B);
104: }
106: /* Test MatTransposeMatMult() */
107: if (!Aiselemental) {
108: MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A^T*A */
109: MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&D);
110: MatTransposeMatMultEqual(A,A,D,10,&equal);
113: /* Test MatDuplicate for matrix product */
114: MatDuplicate(D,MAT_COPY_VALUES,&C);
115: MatDestroy(&C);
116: MatDestroy(&D);
118: /* Test D*x = A^T*C*A*x, where C is in AIJ format */
119: MatGetLocalSize(A,&am,&an);
120: MatCreate(PETSC_COMM_WORLD,&C);
121: if (size == 1) {
122: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,am,am);
123: } else {
124: MatSetSizes(C,am,am,PETSC_DECIDE,PETSC_DECIDE);
125: }
126: MatSetFromOptions(C);
127: MatSetUp(C);
128: MatGetOwnershipRange(C,&rstart,&rend);
129: v[0] = 1.0;
130: for (i=rstart; i<rend; i++) {
131: MatSetValues(C,1,&i,1,&i,v,INSERT_VALUES);
132: }
133: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
134: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
136: /* B = C*A, D = A^T*B */
137: MatMatMult(C,A,MAT_INITIAL_MATRIX,1.0,&B);
138: MatTransposeMatMult(A,B,MAT_INITIAL_MATRIX,fill,&D);
139: MatTransposeMatMultEqual(A,B,D,10,&equal);
142: MatDestroy(&D);
143: MatDestroy(&C);
144: MatDestroy(&B);
145: }
147: /* Test MatMatTransposeMult() */
148: if (!Aiselemental) {
149: PetscReal diff, scale;
150: PetscInt am, an, aM, aN;
152: MatGetLocalSize(A, &am, &an);
153: MatGetSize(A, &aM, &aN);
154: MatCreateDense(PetscObjectComm((PetscObject)A),PETSC_DECIDE, an, aM + 10, aN, NULL, &B);
155: MatSetRandom(B, NULL);
156: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
157: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
158: MatMatTransposeMult(A,B,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
160: /* Test MatDuplicate for matrix product */
161: MatDuplicate(D,MAT_COPY_VALUES,&C);
163: MatMatTransposeMult(A,B,MAT_REUSE_MATRIX,fill,&D);
164: MatAXPY(C, -1., D, SAME_NONZERO_PATTERN);
166: MatNorm(C, NORM_FROBENIUS, &diff);
167: MatNorm(D, NORM_FROBENIUS, &scale);
169: MatDestroy(&C);
171: MatMatTransposeMultEqual(A,B,D,10,&equal);
173: MatDestroy(&D);
174: MatDestroy(&B);
176: }
178: MatDestroy(&A);
179: PetscFree(v);
180: PetscFinalize();
181: return 0;
182: }
184: /*TEST
186: test:
187: output_file: output/ex104.out
189: test:
190: suffix: 2
191: nsize: 2
192: output_file: output/ex104.out
194: test:
195: suffix: 3
196: nsize: 4
197: output_file: output/ex104.out
198: args: -M 23 -N 31
200: test:
201: suffix: 4
202: nsize: 4
203: output_file: output/ex104.out
204: args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via cyclic
206: test:
207: suffix: 5
208: nsize: 4
209: output_file: output/ex104.out
210: args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via allgatherv
212: test:
213: suffix: 6
214: args: -mat_type elemental
215: requires: elemental
216: output_file: output/ex104.out
218: test:
219: suffix: 7
220: nsize: 2
221: args: -mat_type elemental
222: requires: elemental
223: output_file: output/ex104.out
225: TEST*/