Actual source code: ex54.c

  1: /*$Id: ex54.c,v 1.22 2001/08/07 21:30:08 bsmith Exp $*/

  3: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for MatBAIJ format.\n";

 5:  #include petscmat.h

  9: int main(int argc,char **args)
 10: {
 11:   Mat         A,B,*submatA,*submatB;
 12:   int         bs=1,m=11,ov=1,i,j,k,*rows,*cols,ierr,nd=5,*idx,size;
 13:   int         rank,rstart,rend,sz,mm,nn,M,N,Mbs;
 14:   PetscScalar *vals,rval;
 15:   IS          *is1,*is2;
 16:   PetscRandom rand;
 17:   Vec         xx,s1,s2;
 18:   PetscReal   s1norm,s2norm,rnorm,tol = 1.e-10;
 19:   PetscTruth  flg;

 21:   PetscInitialize(&argc,&args,(char *)0,help);
 22:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 23:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 25:   PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
 26:   PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);
 27:   PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
 28:   PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);

 30:   MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,
 31:                           PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL,&A);
 32:   MatCreateMPIAIJ(PETSC_COMM_WORLD,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,
 33:                          PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL,&B);
 34:   PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rand);

 36:   MatGetOwnershipRange(A,&rstart,&rend);
 37:   MatGetSize(A,&M,&N);
 38:   Mbs  = M/bs;

 40:   PetscMalloc(bs*sizeof(int),&rows);
 41:   PetscMalloc(bs*sizeof(int),&cols);
 42:   PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
 43:   PetscMalloc(M*sizeof(PetscScalar),&idx);

 45:   /* Now set blocks of values */
 46:   for (i=0; i<40*bs; i++) {
 47:       PetscRandomGetValue(rand,&rval);
 48:       cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
 49:       PetscRandomGetValue(rand,&rval);
 50:       rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m);
 51:       for (j=1; j<bs; j++) {
 52:         rows[j] = rows[j-1]+1;
 53:         cols[j] = cols[j-1]+1;
 54:       }

 56:       for (j=0; j<bs*bs; j++) {
 57:         PetscRandomGetValue(rand,&rval);
 58:         vals[j] = rval;
 59:       }
 60:       MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
 61:       MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
 62:   }

 64:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 65:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 67:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 69:     /* Test MatIncreaseOverlap() */
 70:   PetscMalloc(nd*sizeof(IS **),&is1);
 71:   PetscMalloc(nd*sizeof(IS **),&is2);

 73: 
 74:   for (i=0; i<nd; i++) {
 75:     PetscRandomGetValue(rand,&rval);
 76:     sz = (int)(PetscRealPart(rval)*m);
 77:     for (j=0; j<sz; j++) {
 78:       PetscRandomGetValue(rand,&rval);
 79:       idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
 80:       for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
 81:     }
 82:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is1+i);
 83:     ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is2+i);
 84:   }
 85:   MatIncreaseOverlap(A,nd,is1,ov);
 86:   MatIncreaseOverlap(B,nd,is2,ov);

 88:   for (i=0; i<nd; ++i) {
 89:     ISEqual(is1[i],is2[i],&flg);

 91:     if (!flg) {
 92:       PetscPrintf(PETSC_COMM_SELF,"i=%d, flg=%d :bs=%d m=%d ov=%d nd=%d np=%d\n",i,flg,bs,m,ov,nd,size);
 93:     }
 94:   }

 96:   for (i=0; i<nd; ++i) {
 97:     ISSort(is1[i]);
 98:     ISSort(is2[i]);
 99:   }
100: 
101:   MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
102:   MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);


105:   /* Test MatMult() */
106:   for (i=0; i<nd; i++) {
107:     MatGetSize(submatA[i],&mm,&nn);
108:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
109:     VecDuplicate(xx,&s1);
110:     VecDuplicate(xx,&s2);
111:     for (j=0; j<3; j++) {
112:       VecSetRandom(rand,xx);
113:       MatMult(submatA[i],xx,s1);
114:       MatMult(submatB[i],xx,s2);
115:       VecNorm(s1,NORM_2,&s1norm);
116:       VecNorm(s2,NORM_2,&s2norm);
117:       rnorm = s2norm-s1norm;
118:       if (rnorm<-tol || rnorm>tol) {
119:         PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
120:       }
121:     }
122:     VecDestroy(xx);
123:     VecDestroy(s1);
124:     VecDestroy(s2);
125:   }

127:   /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
128: 
129:   MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
130:   MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);

132:   /* Test MatMult() */
133:   for (i=0; i<nd; i++) {
134:     MatGetSize(submatA[i],&mm,&nn);
135:     VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
136:     VecDuplicate(xx,&s1);
137:     VecDuplicate(xx,&s2);
138:     for (j=0; j<3; j++) {
139:       VecSetRandom(rand,xx);
140:       MatMult(submatA[i],xx,s1);
141:       MatMult(submatB[i],xx,s2);
142:       VecNorm(s1,NORM_2,&s1norm);
143:       VecNorm(s2,NORM_2,&s2norm);
144:       rnorm = s2norm-s1norm;
145:       if (rnorm<-tol || rnorm>tol) {
146:         PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
147:       }
148:     }
149:     VecDestroy(xx);
150:     VecDestroy(s1);
151:     VecDestroy(s2);
152:   }
153: 
154:   /* Free allocated memory */
155:   for (i=0; i<nd; ++i) {
156:     ISDestroy(is1[i]);
157:     ISDestroy(is2[i]);
158:     MatDestroy(submatA[i]);
159:     MatDestroy(submatB[i]);
160:  }
161:   PetscFree(is1);
162:   PetscFree(is2);
163:   PetscFree(idx);
164:   PetscFree(rows);
165:   PetscFree(cols);
166:   PetscFree(vals);
167:   MatDestroy(A);
168:   MatDestroy(B);
169:   PetscFree(submatA);
170:   PetscFree(submatB);
171:   PetscRandomDestroy(rand);
172:   PetscFinalize();
173:   return 0;
174: }