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;
 17:   PetscReal      nrm;
 19:   PetscInt       size=8,size1=6,size2=2, i,j;

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

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

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

 48:   PetscMalloc(size*sizeof(PetscScalar),&y);
 49:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,y,&Y);
 50:   VecAssemblyBegin(Y);
 51:   VecAssemblyEnd(Y);

 53:   PetscMalloc(size*sizeof(PetscScalar),&z);
 54:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,z,&Z);
 55:   VecAssemblyBegin(Z);
 56:   VecAssemblyEnd(Z);

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

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

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

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

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

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

111:   PetscPrintf(PETSC_COMM_WORLD,"MatMult the usual way:\n");
112:   VecView(Y,0);
113:   PetscPrintf(PETSC_COMM_WORLD,"MatMult by subblock:\n");
114:   VecView(Z,0);

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

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

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

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

166:   VecDestroy(&X);
167:   VecDestroy(&Y);
168:   VecDestroy(&Z);

170:   VecDestroy(&X1);
171:   VecDestroy(&X2);
172:   VecDestroy(&Z1);
173:   VecDestroy(&Z2);

175:   PetscFinalize();
176:   return 0;
177: }