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*/