Actual source code: ex2.c

  1: /*
  2:  * Example code testing SeqDense matrices with an LDA (leading dimension
  3:  * of the user-allocated arrray) larger than M.
  4:  * This example tests the functionality of MatInsertValues, MatMult,
  5:  * and MatMultTranspose.
  6:  */
  7: #include <stdlib.h>
 8:  #include petscmat.h

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

 20:   PetscInitialize(&argc,&argv,0,0);

 22:   /*
 23:    * Create matrix and three vectors: these are all normal
 24:    */
 25:   PetscMalloc(size*size*sizeof(PetscScalar),&a);
 26:   PetscMalloc(size*size*sizeof(PetscScalar),&b);
 27:   for (i=0; i<size; i++) {
 28:     for (j=0; j<size; j++) {
 29:       a[i+j*size] = rand(); b[i+j*size] = a[i+j*size];
 30:     }
 31:   }
 32:   MatCreate(MPI_COMM_SELF,size,size,size,size,&A);
 33:   MatSetType(A,MATSEQDENSE);
 34:   MatSeqDenseSetPreallocation(A,a);
 35:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 36:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 38:   PetscMalloc(size*sizeof(PetscScalar),&x);
 39:   for (i=0; i<size; i++) {
 40:     x[i] = one;
 41:   }
 42:   VecCreateSeqWithArray(MPI_COMM_SELF,size,x,&X);
 43:   VecAssemblyBegin(X);
 44:   VecAssemblyEnd(X);

 46:   PetscMalloc(size*sizeof(PetscScalar),&y);
 47:   VecCreateSeqWithArray(MPI_COMM_SELF,size,y,&Y);
 48:   VecAssemblyBegin(Y);
 49:   VecAssemblyEnd(Y);

 51:   PetscMalloc(size*sizeof(PetscScalar),&z);
 52:   VecCreateSeqWithArray(MPI_COMM_SELF,size,z,&Z);
 53:   VecAssemblyBegin(Z);
 54:   VecAssemblyEnd(Z);

 56:   /*
 57:    * Now create submatrices and subvectors
 58:    */
 59:   MatCreate(MPI_COMM_SELF,size1,size1,size1,size1,&A11);
 60:   MatSetType(A11,MATSEQDENSE);
 61:   MatSeqDenseSetPreallocation(A11,b);
 62:   MatSeqDenseSetLDA(A11,size);
 63:   MatAssemblyBegin(A11,MAT_FINAL_ASSEMBLY);
 64:   MatAssemblyEnd(A11,MAT_FINAL_ASSEMBLY);

 66:   MatCreate(MPI_COMM_SELF,size1,size2,size1,size2,&A12);
 67:   MatSetType(A12,MATSEQDENSE);
 68:   MatSeqDenseSetPreallocation(A12,b+size1*size);
 69:   MatSeqDenseSetLDA(A12,size);
 70:   MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
 71:   MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);

 73:   MatCreate(MPI_COMM_SELF,size2,size1,size2,size1,&A21);
 74:   MatSetType(A21,MATSEQDENSE);
 75:   MatSeqDenseSetPreallocation(A21,b+size1);
 76:   MatSeqDenseSetLDA(A21,size);
 77:   MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);
 78:   MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);

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

 87:   VecCreateSeqWithArray(MPI_COMM_SELF,size1,x,&X1);
 88:   VecCreateSeqWithArray(MPI_COMM_SELF,size2,x+size1,&X2);
 89:   VecCreateSeqWithArray(MPI_COMM_SELF,size1,z,&Z1);
 90:   VecCreateSeqWithArray(MPI_COMM_SELF,size2,z+size1,&Z2);

 92:   /*
 93:    * Now multiple matrix times input in two ways;
 94:    * compare the results
 95:    */
 96:   MatMult(A,X,Y);
 97:   MatMult(A11,X1,Z1);
 98:   MatMultAdd(A12,X2,Z1,Z1);
 99:   MatMult(A22,X2,Z2);
100:   MatMultAdd(A21,X1,Z2,Z2);
101:   VecAXPY(&mone,Y,Z);
102:   VecNorm(Z,NORM_2,&nrm);
103:   printf("Test1; error norm=%e\n",nrm);
104:   /*
105:   printf("MatMult the usual way:\n"); VecView(Y,0);
106:   printf("MatMult by subblock:\n"); VecView(Z,0);
107:   */

109:   /*
110:    * Next test: change both matrices
111:    */
112:   v = rand(); i=1; j=size-2;
113:   MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
114:   j -= size1;
115:   MatSetValues(A12,1,&i,1,&j,&v,INSERT_VALUES);
116:   v = rand(); i=j=size1+1;
117:   MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
118:   i=j=1;
119:   MatSetValues(A22,1,&i,1,&j,&v,INSERT_VALUES);
120:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
121:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
122:   MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);
123:   MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);
124:   MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
125:   MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);

127:   MatMult(A,X,Y);
128:   MatMult(A11,X1,Z1);
129:   MatMultAdd(A12,X2,Z1,Z1);
130:   MatMult(A22,X2,Z2);
131:   MatMultAdd(A21,X1,Z2,Z2);
132:   VecAXPY(&mone,Y,Z);
133:   VecNorm(Z,NORM_2,&nrm);
134:   printf("Test2; error norm=%e\n",nrm);

136:   /*
137:    * Transpose product
138:    */
139:   MatMultTranspose(A,X,Y);
140:   MatMultTranspose(A11,X1,Z1);
141:   MatMultTransposeAdd(A21,X2,Z1,Z1);
142:   MatMultTranspose(A22,X2,Z2);
143:   MatMultTransposeAdd(A12,X1,Z2,Z2);
144:   VecAXPY(&mone,Y,Z);
145:   VecNorm(Z,NORM_2,&nrm);
146:   printf("Test3; error norm=%e\n",nrm);

148:   PetscFinalize();
149:   return 0;
150: }