Actual source code: ex104.c
petsc-3.13.6 2020-09-29
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;
14: PetscRandom r;
15: PetscBool equal,Aiselemental;
16: PetscReal fill = 1.0;
17: IS isrows,iscols;
18: const PetscInt *rows,*cols;
19: PetscScalar *v,rval;
20: #if defined(PETSC_HAVE_ELEMENTAL)
21: PetscBool Test_MatMatMult=PETSC_TRUE;
22: #else
23: PetscBool Test_MatMatMult=PETSC_FALSE;
24: #endif
25: PetscMPIInt size;
27: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
28: MPI_Comm_size(PETSC_COMM_WORLD,&size);
30: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
31: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
32: MatCreate(PETSC_COMM_WORLD,&A);
33: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
34: MatSetType(A,MATDENSE);
35: MatSetFromOptions(A);
36: MatSetUp(A);
37: PetscRandomCreate(PETSC_COMM_WORLD,&r);
38: PetscRandomSetFromOptions(r);
40: /* Set local matrix entries */
41: MatGetOwnershipIS(A,&isrows,&iscols);
42: ISGetLocalSize(isrows,&nrows);
43: ISGetIndices(isrows,&rows);
44: ISGetLocalSize(iscols,&ncols);
45: ISGetIndices(iscols,&cols);
46: PetscMalloc1(nrows*ncols,&v);
47: for (i=0; i<nrows; i++) {
48: for (j=0; j<ncols; j++) {
49: PetscRandomGetValue(r,&rval);
50: v[i*ncols+j] = rval;
51: }
52: }
53: MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES);
54: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
55: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
56: ISRestoreIndices(isrows,&rows);
57: ISRestoreIndices(iscols,&cols);
58: ISDestroy(&isrows);
59: ISDestroy(&iscols);
60: PetscRandomDestroy(&r);
62: PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&Aiselemental);
64: /* Test MatCreateTranspose() and MatTranspose() */
65: MatCreateTranspose(A,&C);
66: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
67: MatMultEqual(C,B,10,&equal);
68: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
69: MatDestroy(&B);
71: MatDuplicate(A,MAT_COPY_VALUES,&B);
72: if (!Aiselemental) {
73: MatTranspose(B,MAT_INPLACE_MATRIX,&B);
74: MatMultEqual(C,B,10,&equal);
75: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"C*x != B*x");
76: }
77: MatDestroy(&B);
79: /* Test B = C*A for matrix type transpose and seqdense */
80: if (size == 1 && !Aiselemental) {
81: MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&B);
82: MatMatMultEqual(C,A,B,10,&equal);
83: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B != C*A for matrix type transpose and seqdense");
84: MatDestroy(&B);
85: }
86: MatDestroy(&C);
88: /* Test MatMatMult() */
89: if (Test_MatMatMult) {
90: #if !defined(PETSC_HAVE_ELEMENTAL)
91: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This test requires ELEMENTAL");
92: #endif
93: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
94: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */
95: MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C);
96: MatMatMultEqual(B,A,C,10,&equal);
97: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
99: /* Test MatDuplicate for matrix product */
100: MatDuplicate(C,MAT_COPY_VALUES,&D);
102: MatDestroy(&D);
103: MatDestroy(&C);
104: MatDestroy(&B);
105: }
107: /* Test MatTransposeMatMult() */
108: if (!Aiselemental) {
109: MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A^T*A */
110: MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&D);
111: MatTransposeMatMultEqual(A,A,D,10,&equal);
112: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x");
114: /* Test MatDuplicate for matrix product */
115: MatDuplicate(D,MAT_COPY_VALUES,&C);
116: MatDestroy(&C);
117: MatDestroy(&D);
119: /* Test D*x = A^T*C*A*x, where C is in AIJ format */
120: MatGetLocalSize(A,&am,&an);
121: MatCreate(PETSC_COMM_WORLD,&C);
122: if (size == 1) {
123: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,am,am);
124: } else {
125: MatSetSizes(C,am,am,PETSC_DECIDE,PETSC_DECIDE);
126: }
127: MatSetFromOptions(C);
128: MatSetUp(C);
129: MatGetOwnershipRange(C,&rstart,&rend);
130: v[0] = 1.0;
131: for (i=rstart; i<rend; i++) {
132: MatSetValues(C,1,&i,1,&i,v,INSERT_VALUES);
133: }
134: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
135: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
137: /* B = C*A, D = A^T*B */
138: MatMatMult(C,A,MAT_INITIAL_MATRIX,1.0,&B);
139: MatTransposeMatMult(A,B,MAT_INITIAL_MATRIX,fill,&D);
140: MatTransposeMatMultEqual(A,B,D,10,&equal);
141: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*B*x");
143: MatDestroy(&D);
144: MatDestroy(&C);
145: MatDestroy(&B);
146: }
148: /* Test MatMatTransposeMult() */
149: if (!Aiselemental) {
150: PetscReal diff, scale;
151: PetscInt am, an, aM, aN;
153: MatGetLocalSize(A, &am, &an);
154: MatGetSize(A, &aM, &aN);
155: MatCreateDense(PetscObjectComm((PetscObject)A),PETSC_DECIDE, an, aM + 10, aN, NULL, &B);
156: MatSetRandom(B, NULL);
157: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
158: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
159: MatMatTransposeMult(A,B,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
161: /* Test MatDuplicate for matrix product */
162: MatDuplicate(D,MAT_COPY_VALUES,&C);
164: MatMatTransposeMult(A,B,MAT_REUSE_MATRIX,fill,&D);
165: MatAXPY(C, -1., D, SAME_NONZERO_PATTERN);
167: MatNorm(C, NORM_FROBENIUS, &diff);
168: MatNorm(D, NORM_FROBENIUS, &scale);
169: if (diff > PETSC_SMALL * scale) SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatMatTransposeMult() differs between MAT_INITIAL_MATRIX and MAT_REUSE_MATRIX");
170: MatDestroy(&C);
172: MatMatTransposeMultEqual(A,B,D,10,&equal);
173: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x");
174: MatDestroy(&D);
175: MatDestroy(&B);
177: }
179: MatDestroy(&A);
180: PetscFree(v);
181: PetscFinalize();
182: return ierr;
183: }
185: /*TEST
187: test:
188: output_file: output/ex104.out
190: test:
191: suffix: 2
192: nsize: 2
193: output_file: output/ex104.out
195: test:
196: suffix: 3
197: nsize: 4
198: output_file: output/ex104.out
199: args: -M 23 -N 31
201: test:
202: suffix: 4
203: nsize: 4
204: output_file: output/ex104.out
205: args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via cyclic
207: test:
208: suffix: 5
209: nsize: 4
210: output_file: output/ex104.out
211: args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via allgatherv
213: test:
214: suffix: 6
215: args: -mat_type elemental
216: requires: elemental
217: output_file: output/ex104.out
219: test:
220: suffix: 7
221: nsize: 2
222: args: -mat_type elemental
223: requires: elemental
224: output_file: output/ex104.out
226: TEST*/