Actual source code: ex99.c
1: static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, square, shifted (copied from ex97)\n";
3: #include <petscmat.h>
5: static PetscErrorCode AssembleMatrix(MPI_Comm comm,Mat *A)
6: {
7: Mat B;
8: PetscInt i,ms,me;
10: MatCreate(comm,&B);
11: MatSetSizes(B,6,6,PETSC_DETERMINE,PETSC_DETERMINE);
12: MatSetFromOptions(B);
13: MatSetUp(B);
14: MatGetOwnershipRange(B,&ms,&me);
15: for (i=ms; i<me; i++) {
16: MatSetValue(B,i,i,1.0*i,INSERT_VALUES);
17: }
18: MatSetValue(B,me-1,me-1,me*me,INSERT_VALUES);
19: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
20: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
21: *A = B;
22: return 0;
23: }
25: static PetscErrorCode Compare2(Vec *X,const char *test)
26: {
27: PetscReal norm;
28: Vec Y;
29: PetscInt verbose = 0;
31: VecDuplicate(X[0],&Y);
32: VecCopy(X[0],Y);
33: VecAYPX(Y,-1.0,X[1]);
34: VecNorm(Y,NORM_INFINITY,&norm);
36: PetscOptionsGetInt(NULL,NULL,"-verbose",&verbose,NULL);
37: if (norm < PETSC_SQRT_MACHINE_EPSILON && verbose < 1) {
38: PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference < sqrt(eps_machine)\n",test);
39: } else {
40: PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference %g\n",test,(double)norm);
41: }
42: if (verbose > 1) {
43: VecView(X[0],PETSC_VIEWER_STDOUT_WORLD);
44: VecView(X[1],PETSC_VIEWER_STDOUT_WORLD);
45: VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
46: }
47: VecDestroy(&Y);
48: return 0;
49: }
51: static PetscErrorCode CheckMatrices(Mat A,Mat B,Vec left,Vec right,Vec X,Vec Y,Vec X1,Vec Y1)
52: {
53: Vec *ltmp,*rtmp;
55: VecDuplicateVecs(right,2,&rtmp);
56: VecDuplicateVecs(left,2,<mp);
57: MatScale(A,PETSC_PI);
58: MatScale(B,PETSC_PI);
59: MatDiagonalScale(A,left,right);
60: MatDiagonalScale(B,left,right);
61: MatShift(A,PETSC_PI);
62: MatShift(B,PETSC_PI);
64: MatMult(A,X,ltmp[0]);
65: MatMult(B,X,ltmp[1]);
66: Compare2(ltmp,"MatMult");
68: MatMultTranspose(A,Y,rtmp[0]);
69: MatMultTranspose(B,Y,rtmp[1]);
70: Compare2(rtmp,"MatMultTranspose");
72: VecCopy(Y1,ltmp[0]);
73: VecCopy(Y1,ltmp[1]);
74: MatMultAdd(A,X,ltmp[0],ltmp[0]);
75: MatMultAdd(B,X,ltmp[1],ltmp[1]);
76: Compare2(ltmp,"MatMultAdd v2==v3");
78: MatMultAdd(A,X,Y1,ltmp[0]);
79: MatMultAdd(B,X,Y1,ltmp[1]);
80: Compare2(ltmp,"MatMultAdd v2!=v3");
82: VecCopy(X1,rtmp[0]);
83: VecCopy(X1,rtmp[1]);
84: MatMultTransposeAdd(A,Y,rtmp[0],rtmp[0]);
85: MatMultTransposeAdd(B,Y,rtmp[1],rtmp[1]);
86: Compare2(rtmp,"MatMultTransposeAdd v2==v3");
88: MatMultTransposeAdd(A,Y,X1,rtmp[0]);
89: MatMultTransposeAdd(B,Y,X1,rtmp[1]);
90: Compare2(rtmp,"MatMultTransposeAdd v2!=v3");
92: VecDestroyVecs(2,<mp);
93: VecDestroyVecs(2,&rtmp);
94: return 0;
95: }
97: int main(int argc, char *argv[])
98: {
99: Mat A,B,Asub,Bsub;
100: PetscInt ms,idxrow[3],idxcol[3];
101: Vec left,right,X,Y,X1,Y1;
102: IS isrow,iscol;
103: PetscBool random = PETSC_TRUE;
105: PetscInitialize(&argc,&argv,NULL,help);
106: AssembleMatrix(PETSC_COMM_WORLD,&A);
107: AssembleMatrix(PETSC_COMM_WORLD,&B);
108: MatSetOperation(B,MATOP_CREATE_SUBMATRIX,NULL);
109: MatSetOperation(B,MATOP_CREATE_SUBMATRICES,NULL);
110: MatGetOwnershipRange(A,&ms,NULL);
112: idxrow[0] = ms+1;
113: idxrow[1] = ms+2;
114: idxrow[2] = ms+4;
115: ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow);
117: idxcol[0] = ms+1;
118: idxcol[1] = ms+2;
119: idxcol[2] = ms+4;
120: ISCreateGeneral(PETSC_COMM_WORLD,3,idxcol,PETSC_USE_POINTER,&iscol);
122: MatCreateSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub);
123: MatCreateSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub);
125: MatCreateVecs(Asub,&right,&left);
126: VecDuplicate(right,&X);
127: VecDuplicate(right,&X1);
128: VecDuplicate(left,&Y);
129: VecDuplicate(left,&Y1);
131: PetscOptionsGetBool(NULL,NULL,"-random",&random,NULL);
132: if (random) {
133: VecSetRandom(right,NULL);
134: VecSetRandom(left,NULL);
135: VecSetRandom(X,NULL);
136: VecSetRandom(Y,NULL);
137: VecSetRandom(X1,NULL);
138: VecSetRandom(Y1,NULL);
139: } else {
140: VecSet(right,1.0);
141: VecSet(left,2.0);
142: VecSet(X,3.0);
143: VecSet(Y,4.0);
144: VecSet(X1,3.0);
145: VecSet(Y1,4.0);
146: }
147: CheckMatrices(Asub,Bsub,left,right,X,Y,X1,Y1);
148: ISDestroy(&isrow);
149: ISDestroy(&iscol);
150: MatDestroy(&A);
151: MatDestroy(&B);
152: MatDestroy(&Asub);
153: MatDestroy(&Bsub);
154: VecDestroy(&left);
155: VecDestroy(&right);
156: VecDestroy(&X);
157: VecDestroy(&Y);
158: VecDestroy(&X1);
159: VecDestroy(&Y1);
160: PetscFinalize();
161: return 0;
162: }
164: /*TEST
166: test:
167: nsize: 3
169: TEST*/