Actual source code: ex2.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1: static char help[] = "testing SeqDense matrices with an LDA (leading dimension of the user-allocated arrray) larger than M.\n";
  2: /*
  3:  * Example code testing SeqDense matrices with an LDA (leading dimension
  4:  * of the user-allocated arrray) larger than M.
  5:  * This example tests the functionality of MatSetValues(), MatMult(),
  6:  * and MatMultTranspose().
  7:  */

  9:  #include <petscmat.h>

 11: int main(int argc,char **argv)
 12: {
 13:   Mat            A,A11,A12,A21,A22;
 14:   Vec            X,X1,X2,Y,Z,Z1,Z2;
 15:   PetscScalar    *a,*b,*x,*y,*z,v,one=1;
 16:   PetscReal      nrm;
 18:   PetscInt       size=8,size1=6,size2=2, i,j;
 19:   PetscRandom    rnd;

 21:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
 22:   PetscRandomCreate(PETSC_COMM_SELF,&rnd);

 24:   /*
 25:    * Create matrix and three vectors: these are all normal
 26:    */
 27:   PetscMalloc1(size*size,&a);
 28:   PetscMalloc1(size*size,&b);
 29:   for (i=0; i<size; i++) {
 30:     for (j=0; j<size; j++) {
 31:       PetscRandomGetValue(rnd,&a[i+j*size]);
 32:       b[i+j*size] = a[i+j*size];
 33:     }
 34:   }
 35:   MatCreate(MPI_COMM_SELF,&A);
 36:   MatSetSizes(A,size,size,size,size);
 37:   MatSetType(A,MATSEQDENSE);
 38:   MatSeqDenseSetPreallocation(A,a);
 39:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 40:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 42:   PetscMalloc1(size,&x);
 43:   for (i=0; i<size; i++) {
 44:     x[i] = one;
 45:   }
 46:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,x,&X);
 47:   VecAssemblyBegin(X);
 48:   VecAssemblyEnd(X);

 50:   PetscMalloc1(size,&y);
 51:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,y,&Y);
 52:   VecAssemblyBegin(Y);
 53:   VecAssemblyEnd(Y);

 55:   PetscMalloc1(size,&z);
 56:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,z,&Z);
 57:   VecAssemblyBegin(Z);
 58:   VecAssemblyEnd(Z);

 60:   /*
 61:    * Now create submatrices and subvectors
 62:    */
 63:   MatCreate(MPI_COMM_SELF,&A11);
 64:   MatSetSizes(A11,size1,size1,size1,size1);
 65:   MatSetType(A11,MATSEQDENSE);
 66:   MatSeqDenseSetPreallocation(A11,b);
 67:   MatSeqDenseSetLDA(A11,size);
 68:   MatAssemblyBegin(A11,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(A11,MAT_FINAL_ASSEMBLY);

 71:   MatCreate(MPI_COMM_SELF,&A12);
 72:   MatSetSizes(A12,size1,size2,size1,size2);
 73:   MatSetType(A12,MATSEQDENSE);
 74:   MatSeqDenseSetPreallocation(A12,b+size1*size);
 75:   MatSeqDenseSetLDA(A12,size);
 76:   MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
 77:   MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);

 79:   MatCreate(MPI_COMM_SELF,&A21);
 80:   MatSetSizes(A21,size2,size1,size2,size1);
 81:   MatSetType(A21,MATSEQDENSE);
 82:   MatSeqDenseSetPreallocation(A21,b+size1);
 83:   MatSeqDenseSetLDA(A21,size);
 84:   MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);
 85:   MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);

 87:   MatCreate(MPI_COMM_SELF,&A22);
 88:   MatSetSizes(A22,size2,size2,size2,size2);
 89:   MatSetType(A22,MATSEQDENSE);
 90:   MatSeqDenseSetPreallocation(A22,b+size1*size+size1);
 91:   MatSeqDenseSetLDA(A22,size);
 92:   MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
 93:   MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);

 95:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size1,x,&X1);
 96:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size2,x+size1,&X2);
 97:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size1,z,&Z1);
 98:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size2,z+size1,&Z2);

100:   /*
101:    * Now multiple matrix times input in two ways;
102:    * compare the results
103:    */
104:   MatMult(A,X,Y);
105:   MatMult(A11,X1,Z1);
106:   MatMultAdd(A12,X2,Z1,Z1);
107:   MatMult(A22,X2,Z2);
108:   MatMultAdd(A21,X1,Z2,Z2);
109:   VecAXPY(Z,-1.0,Y);
110:   VecNorm(Z,NORM_2,&nrm);
111:   if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
112:     PetscPrintf(PETSC_COMM_WORLD,"Test1; error norm=%g\n",(double)nrm);
113:   }

115:   /*
116:    * Next test: change both matrices
117:    */
118:   PetscRandomGetValue(rnd,&v);
119:   i    = 1;
120:   j    = size-2;
121:   MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
122:   j   -= size1;
123:   MatSetValues(A12,1,&i,1,&j,&v,INSERT_VALUES);
124:   PetscRandomGetValue(rnd,&v);
125:   i    = j = size1+1;
126:   MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
127:   i     =j = 1;
128:   MatSetValues(A22,1,&i,1,&j,&v,INSERT_VALUES);
129:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
130:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
131:   MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
132:   MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
133:   MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
134:   MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);

136:   MatMult(A,X,Y);
137:   MatMult(A11,X1,Z1);
138:   MatMultAdd(A12,X2,Z1,Z1);
139:   MatMult(A22,X2,Z2);
140:   MatMultAdd(A21,X1,Z2,Z2);
141:   VecAXPY(Z,-1.0,Y);
142:   VecNorm(Z,NORM_2,&nrm);
143:   if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
144:     PetscPrintf(PETSC_COMM_WORLD,"Test2; error norm=%g\n",(double)nrm);
145:   }

147:   /*
148:    * Transpose product
149:    */
150:   MatMultTranspose(A,X,Y);
151:   MatMultTranspose(A11,X1,Z1);
152:   MatMultTransposeAdd(A21,X2,Z1,Z1);
153:   MatMultTranspose(A22,X2,Z2);
154:   MatMultTransposeAdd(A12,X1,Z2,Z2);
155:   VecAXPY(Z,-1.0,Y);
156:   VecNorm(Z,NORM_2,&nrm);
157:   if (nrm > 100.0*PETSC_MACHINE_EPSILON) {
158:     PetscPrintf(PETSC_COMM_WORLD,"Test3; error norm=%g\n",(double)nrm);
159:   }
160:   PetscFree(a);
161:   PetscFree(b);
162:   PetscFree(x);
163:   PetscFree(y);
164:   PetscFree(z);
165:   MatDestroy(&A);
166:   MatDestroy(&A11);
167:   MatDestroy(&A12);
168:   MatDestroy(&A21);
169:   MatDestroy(&A22);

171:   VecDestroy(&X);
172:   VecDestroy(&Y);
173:   VecDestroy(&Z);

175:   VecDestroy(&X1);
176:   VecDestroy(&X2);
177:   VecDestroy(&Z1);
178:   VecDestroy(&Z2);
179:   PetscRandomDestroy(&rnd);
180:   PetscFinalize();
181:   return ierr;
182: }


185: /*TEST

187:    test:

189: TEST*/