Actual source code: ex104.c
petsc-3.12.0 2019-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,iselemental;
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: /* Test MatTranspose() */
63: MatCreateTranspose(A,&C);
64: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
65: MatMultEqual(C,B,10,&equal);
66: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
67: MatTranspose(A,MAT_REUSE_MATRIX,&B); /* B = A^T */
68: MatMultEqual(C,B,10,&equal);
69: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
70: MatDestroy(&B);
71: MatDuplicate(A,MAT_COPY_VALUES,&B);
72: MatTranspose(B,MAT_INPLACE_MATRIX,&B);
73: MatMultEqual(C,B,10,&equal);
74: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
75: MatDestroy(&B);
76: MatDestroy(&C);
78: /* Test MatMatMult() */
79: if (Test_MatMatMult) {
80: #if !defined(PETSC_HAVE_ELEMENTAL)
81: if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This test requires ELEMENTAL");
82: #endif
83: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
84: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */
85: MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C);
87: /* Test MatDuplicate for matrix product */
88: MatDuplicate(C,MAT_COPY_VALUES,&D);
89: MatDestroy(&D);
91: /* Test B*A*x = C*x for n random vector x */
92: MatMatMultEqual(B,A,C,10,&equal);
93: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
94: MatDestroy(&C);
96: MatMatMultSymbolic(B,A,fill,&C);
97: for (i=0; i<2; i++) {
98: /* Repeat the numeric product to test reuse of the previous symbolic product */
99: MatMatMultNumeric(B,A,C);
101: MatMatMultEqual(B,A,C,10,&equal);
102: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
103: }
104: MatDestroy(&C);
105: MatDestroy(&B);
106: }
108: /* Test MatTransposeMatMult() */
109: PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&iselemental);
110: if (!iselemental) {
111: MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A^T*A */
112: MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&D);
114: /* Test MatDuplicate for matrix product */
115: MatDuplicate(D,MAT_COPY_VALUES,&C);
116: MatDestroy(&C);
118: MatTransposeMatMultEqual(A,A,D,10,&equal);
119: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x");
120: MatDestroy(&D);
122: /* Test D*x = A^T*C*A*x, where C is in AIJ format */
123: MatGetLocalSize(A,&am,&an);
124: MatCreate(PETSC_COMM_WORLD,&C);
125: if (size == 1) {
126: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,am,am);
127: } else {
128: MatSetSizes(C,am,am,PETSC_DECIDE,PETSC_DECIDE);
129: }
130: MatSetFromOptions(C);
131: MatSetUp(C);
132: MatGetOwnershipRange(C,&rstart,&rend);
133: v[0] = 1.0;
134: for (i=rstart; i<rend; i++) {
135: MatSetValues(C,1,&i,1,&i,v,INSERT_VALUES);
136: }
137: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
138: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
140: /* B = C*A, D = A^T*B */
141: MatMatMult(C,A,MAT_INITIAL_MATRIX,1.0,&B);
142: MatTransposeMatMult(A,B,MAT_INITIAL_MATRIX,fill,&D);
143: MatTransposeMatMultEqual(A,B,D,10,&equal);
144: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*B*x");
146: MatDestroy(&D);
147: MatDestroy(&C);
148: MatDestroy(&B);
149: }
151: /* Test MatMatTransposeMult() */
152: if (!iselemental) {
153: PetscReal diff, scale;
154: PetscInt am, an, aM, aN;
156: MatGetLocalSize(A, &am, &an);
157: MatGetSize(A, &aM, &aN);
158: MatCreateDense(PetscObjectComm((PetscObject)A),PETSC_DECIDE, an, aM + 10, aN, NULL, &B);
159: MatSetRandom(B, NULL);
160: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
161: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
162: MatMatTransposeMult(A,B,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
164: /* Test MatDuplicate for matrix product */
165: MatDuplicate(D,MAT_COPY_VALUES,&C);
167: MatMatTransposeMult(A,B,MAT_REUSE_MATRIX,fill,&D);
168: MatAXPY(C, -1., D, SAME_NONZERO_PATTERN);
170: MatNorm(C, NORM_FROBENIUS, &diff);
171: MatNorm(D, NORM_FROBENIUS, &scale);
172: if (diff > PETSC_SMALL * scale) SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatMatTransposeMult() differs between MAT_INITIAL_MATRIX and MAT_REUSE_MATRIX");
173: MatDestroy(&C);
175: MatMatTransposeMultEqual(A,B,D,10,&equal);
176: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x");
177: MatDestroy(&D);
178: MatDestroy(&B);
180: }
182: MatDestroy(&A);
183: PetscFree(v);
184: PetscFinalize();
185: return ierr;
186: }
188: /*TEST
190: test:
191: output_file: output/ex104.out
193: test:
194: suffix: 2
195: nsize: 2
196: output_file: output/ex104.out
198: test:
199: suffix: 3
200: nsize: 4
201: output_file: output/ex104.out
202: args: -M 23 -N 31
204: test:
205: suffix: 4
206: nsize: 4
207: output_file: output/ex104.out
208: args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via cyclic
210: test:
211: suffix: 5
212: nsize: 4
213: output_file: output/ex104.out
214: args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via allgatherv
216: TEST*/