Actual source code: ex177.c


  2: static char help[] = "Tests various routines in MatKAIJ format.\n";

  4: #include <petscmat.h>
  5: #define IMAX 15

  7: int main(int argc,char **args)
  8: {
  9:   Mat            A,B,TA;
 10:   PetscScalar    *S,*T;
 11:   PetscViewer    fd;
 12:   char           file[PETSC_MAX_PATH_LEN];
 13:   PetscInt       m,n,M,N,p=1,q=1,i,j;
 14:   PetscMPIInt    rank,size;
 15:   PetscBool      flg;

 17:   PetscInitialize(&argc,&args,(char*)0,help);
 18:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 19:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 21:   /* Load AIJ matrix A */
 22:   PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
 24:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 25:   MatCreate(PETSC_COMM_WORLD,&A);
 26:   MatLoad(A,fd);
 27:   PetscViewerDestroy(&fd);

 29:   /* Get dof, then create S and T */
 30:   PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
 31:   PetscOptionsGetInt(NULL,NULL,"-q",&q,NULL);
 32:   PetscMalloc2(p*q,&S,p*q,&T);
 33:   for (i=0; i<p*q; i++) S[i] = 0;

 35:   for (i=0; i<p; i++) {
 36:     for (j=0; j<q; j++) {
 37:       /* Set some random non-zero values */
 38:       S[i+p*j] = ((PetscReal) ((i+1)*(j+1))) / ((PetscReal) (p+q));
 39:       T[i+p*j] = ((PetscReal) ((p-i)+j)) / ((PetscReal) (p*q));
 40:     }
 41:   }

 43:   /* Test KAIJ when both S & T are not NULL */

 45:   /* Create KAIJ matrix TA */
 46:   MatCreateKAIJ(A,p,q,S,T,&TA);
 47:   MatGetLocalSize(A,&m,&n);
 48:   MatGetSize(A,&M,&N);

 50:   MatConvert(TA,MATAIJ,MAT_INITIAL_MATRIX,&B);

 52:   /* Test MatKAIJGetScaledIdentity() */
 53:   MatKAIJGetScaledIdentity(TA,&flg);
 55:   /* Test MatMult() */
 56:   MatMultEqual(TA,B,10,&flg);
 58:   /* Test MatMultAdd() */
 59:   MatMultAddEqual(TA,B,10,&flg);

 62:   MatDestroy(&TA);
 63:   MatDestroy(&B);

 65:   /* Test KAIJ when S is NULL */

 67:   /* Create KAIJ matrix TA */
 68:   MatCreateKAIJ(A,p,q,NULL,T,&TA);
 69:   MatGetLocalSize(A,&m,&n);
 70:   MatGetSize(A,&M,&N);

 72:   MatConvert(TA,MATAIJ,MAT_INITIAL_MATRIX,&B);

 74:   /* Test MatKAIJGetScaledIdentity() */
 75:   MatKAIJGetScaledIdentity(TA,&flg);
 77:   /* Test MatMult() */
 78:   MatMultEqual(TA,B,10,&flg);
 80:   /* Test MatMultAdd() */
 81:   MatMultAddEqual(TA,B,10,&flg);

 84:   MatDestroy(&TA);
 85:   MatDestroy(&B);

 87:   /* Test KAIJ when T is NULL */

 89:   /* Create KAIJ matrix TA */
 90:   MatCreateKAIJ(A,p,q,S,NULL,&TA);
 91:   MatGetLocalSize(A,&m,&n);
 92:   MatGetSize(A,&M,&N);

 94:   MatConvert(TA,MATAIJ,MAT_INITIAL_MATRIX,&B);

 96:   /* Test MatKAIJGetScaledIdentity() */
 97:   MatKAIJGetScaledIdentity(TA,&flg);
 99:   /* Test MatMult() */
100:   MatMultEqual(TA,B,10,&flg);
102:   /* Test MatMultAdd() */
103:   MatMultAddEqual(TA,B,10,&flg);

106:   MatDestroy(&TA);
107:   MatDestroy(&B);

109:   /* Test KAIJ when T is is an identity matrix */

111:   if (p == q) {
112:     for (i=0; i<p; i++) {
113:       for (j=0; j<q; j++) {
114:         if (i==j) T[i+j*p] = 1.0;
115:         else      T[i+j*p] = 0.0;
116:       }
117:     }

119:     /* Create KAIJ matrix TA */
120:     MatCreateKAIJ(A,p,q,S,T,&TA);
121:     MatGetLocalSize(A,&m,&n);
122:     MatGetSize(A,&M,&N);

124:     MatConvert(TA,MATAIJ,MAT_INITIAL_MATRIX,&B);

126:     /* Test MatKAIJGetScaledIdentity() */
127:     MatKAIJGetScaledIdentity(TA,&flg);
129:     /* Test MatMult() */
130:     MatMultEqual(TA,B,10,&flg);
132:     /* Test MatMultAdd() */
133:     MatMultAddEqual(TA,B,10,&flg);

136:     MatDestroy(&TA);
137:     MatDestroy(&B);

139:     MatCreateKAIJ(A,p,q,NULL,T,&TA);
140:     /* Test MatKAIJGetScaledIdentity() */
141:     MatKAIJGetScaledIdentity(TA,&flg);
143:     MatDestroy(&TA);

145:     for (i=0; i<p; i++) {
146:       for (j=0; j<q; j++) {
147:         if (i==j) S[i+j*p] = T[i+j*p] = 2.0;
148:         else      S[i+j*p] = T[i+j*p] = 0.0;
149:       }
150:     }
151:     MatCreateKAIJ(A,p,q,S,T,&TA);
152:     /* Test MatKAIJGetScaledIdentity() */
153:     MatKAIJGetScaledIdentity(TA,&flg);
155:     MatDestroy(&TA);
156:   }

158:   /* Done with all tests */

160:   PetscFree2(S,T);
161:   MatDestroy(&A);
162:   PetscFinalize();
163:   return 0;
164: }

166: /*TEST

168:   build:
169:     requires: !complex

171:   test:
172:     requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
173:     output_file: output/ex176.out
174:     nsize: {{1 2 3 4}}
175:     args: -f ${DATAFILESPATH}/matrices/small -p {{2 3 7}} -q {{3 7}} -viewer_binary_skip_info

177: TEST*/