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 MatInsertValues(), MatMult(),
  6:  * and MatMultTranspose().
  7:  */
  8:  #include <petscmat.h>

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

 19:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;

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

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

 46:   PetscMalloc1(size,&y);
 47:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,y,&Y);
 48:   VecAssemblyBegin(Y);
 49:   VecAssemblyEnd(Y);

 51:   PetscMalloc1(size,&z);
 52:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,z,&Z);
 53:   VecAssemblyBegin(Z);
 54:   VecAssemblyEnd(Z);

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

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

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

 83:   MatCreate(MPI_COMM_SELF,&A22);
 84:   MatSetSizes(A22,size2,size2,size2,size2);
 85:   MatSetType(A22,MATSEQDENSE);
 86:   MatSeqDenseSetPreallocation(A22,b+size1*size+size1);
 87:   MatSeqDenseSetLDA(A22,size);
 88:   MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);
 89:   MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);

 91:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size1,x,&X1);
 92:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size2,x+size1,&X2);
 93:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size1,z,&Z1);
 94:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size2,z+size1,&Z2);

 96:   /*
 97:    * Now multiple matrix times input in two ways;
 98:    * compare the results
 99:    */
100:   MatMult(A,X,Y);
101:   MatMult(A11,X1,Z1);
102:   MatMultAdd(A12,X2,Z1,Z1);
103:   MatMult(A22,X2,Z2);
104:   MatMultAdd(A21,X1,Z2,Z2);
105:   VecAXPY(Z,-1.0,Y);
106:   VecNorm(Z,NORM_2,&nrm);
107:   PetscPrintf(PETSC_COMM_WORLD,"Test1; error norm=%g\n",(double)nrm);
108:   PetscPrintf(PETSC_COMM_WORLD,"MatMult the usual way:\n");
109:   VecView(Y,0);
110:   PetscPrintf(PETSC_COMM_WORLD,"MatMult by subblock:\n");
111:   VecView(Z,0);

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

131:   MatMult(A,X,Y);
132:   MatMult(A11,X1,Z1);
133:   MatMultAdd(A12,X2,Z1,Z1);
134:   MatMult(A22,X2,Z2);
135:   MatMultAdd(A21,X1,Z2,Z2);
136:   VecAXPY(Z,-1.0,Y);
137:   VecNorm(Z,NORM_2,&nrm);
138:   PetscPrintf(PETSC_COMM_WORLD,"Test2; error norm=%g\n",(double)nrm);

140:   /*
141:    * Transpose product
142:    */
143:   MatMultTranspose(A,X,Y);
144:   MatMultTranspose(A11,X1,Z1);
145:   MatMultTransposeAdd(A21,X2,Z1,Z1);
146:   MatMultTranspose(A22,X2,Z2);
147:   MatMultTransposeAdd(A12,X1,Z2,Z2);
148:   VecAXPY(Z,-1.0,Y);
149:   VecNorm(Z,NORM_2,&nrm);
150:   PetscPrintf(PETSC_COMM_WORLD,"Test3; error norm=%g\n",(double)nrm);

152:   PetscFree(a);
153:   PetscFree(b);
154:   PetscFree(x);
155:   PetscFree(y);
156:   PetscFree(z);
157:   MatDestroy(&A);
158:   MatDestroy(&A11);
159:   MatDestroy(&A12);
160:   MatDestroy(&A21);
161:   MatDestroy(&A22);

163:   VecDestroy(&X);
164:   VecDestroy(&Y);
165:   VecDestroy(&Z);

167:   VecDestroy(&X1);
168:   VecDestroy(&X2);
169:   VecDestroy(&Z1);
170:   VecDestroy(&Z2);

172:   PetscFinalize();
173:   return ierr;
174: }