Actual source code: ex97.c
1: static const char help[] = "Tests MatGetSubMatrix with MatSubMatrix versus MatAIJ, non-square\n";
3: #include <petscmat.h>
7: static PetscErrorCode AssembleMatrix(MPI_Comm comm,Mat *A)
8: {
10: Mat B;
11: PetscInt i,ms,me;
14: MatCreate(comm,&B);
15: MatSetSizes(B,5,6,PETSC_DETERMINE,PETSC_DETERMINE);
16: MatSetFromOptions(B);
17: MatGetOwnershipRange(B,&ms,&me);
18: for (i=ms; i<me; i++) {
19: MatSetValue(B,i,i,1.0*i,INSERT_VALUES);
20: }
21: MatSetValue(B,me-1,me,me*me,INSERT_VALUES);
22: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
23: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
24: *A = B;
25: return(0);
26: }
30: static PetscErrorCode Compare2(Vec *X,const char *test)
31: {
33: PetscReal norm;
34: Vec Y;
35: PetscInt verbose = 0;
38: VecDuplicate(X[0],&Y);
39: VecCopy(X[0],Y);
40: VecAYPX(Y,-1.0,X[1]);
41: VecNorm(Y,NORM_INFINITY,&norm);
43: PetscOptionsGetInt(PETSC_NULL,"-verbose",&verbose,PETSC_NULL);
44: if (norm < 1.e-12 && verbose < 1) {
45: PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference < 1e-12\n",test);
46: } else {
47: PetscPrintf(PETSC_COMM_WORLD,"%30s: norm difference %G\n",test,norm);
48: }
49: if (verbose > 1) {
50: VecView(X[0],PETSC_VIEWER_STDOUT_WORLD);
51: VecView(X[1],PETSC_VIEWER_STDOUT_WORLD);
52: VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
53: }
54: VecDestroy(&Y);
55: return(0);
56: }
60: static PetscErrorCode CheckMatrices(Mat A,Mat B,Vec left,Vec right,Vec X,Vec Y,Vec X1,Vec Y1)
61: {
63: Vec *ltmp,*rtmp;
66: VecDuplicateVecs(right,2,&rtmp);
67: VecDuplicateVecs(left,2,<mp);
68: MatScale(A,PETSC_PI);
69: MatScale(B,PETSC_PI);
70: MatDiagonalScale(A,left,right);
71: MatDiagonalScale(B,left,right);
73: MatMult(A,X,ltmp[0]);
74: MatMult(B,X,ltmp[1]);
75: Compare2(ltmp,"MatMult");
77: MatMultTranspose(A,Y,rtmp[0]);
78: MatMultTranspose(B,Y,rtmp[1]);
79: Compare2(rtmp,"MatMultTranspose");
81: VecCopy(Y1,ltmp[0]);
82: VecCopy(Y1,ltmp[1]);
83: MatMultAdd(A,X,ltmp[0],ltmp[0]);
84: MatMultAdd(B,X,ltmp[1],ltmp[1]);
85: Compare2(ltmp,"MatMultAdd v2==v3");
87: MatMultAdd(A,X,Y1,ltmp[0]);
88: MatMultAdd(B,X,Y1,ltmp[1]);
89: Compare2(ltmp,"MatMultAdd v2!=v3");
91: VecCopy(X1,rtmp[0]);
92: VecCopy(X1,rtmp[1]);
93: MatMultTransposeAdd(A,Y,rtmp[0],rtmp[0]);
94: MatMultTransposeAdd(B,Y,rtmp[1],rtmp[1]);
95: Compare2(rtmp,"MatMultTransposeAdd v2==v3");
97: MatMultTransposeAdd(A,Y,X1,rtmp[0]);
98: MatMultTransposeAdd(B,Y,X1,rtmp[1]);
99: Compare2(rtmp,"MatMultTransposeAdd v2!=v3");
101: VecDestroyVecs(2,<mp);
102: VecDestroyVecs(2,&rtmp);
103: return(0);
104: }
108: int main(int argc, char *argv[])
109: {
111: Mat A,B,Asub,Bsub;
112: PetscInt ms,idxrow[3],idxcol[4];
113: Vec left,right,X,Y,X1,Y1;
114: IS isrow,iscol;
115: PetscBool random = PETSC_TRUE;
117: PetscInitialize(&argc,&argv,PETSC_NULL,help);
118: AssembleMatrix(PETSC_COMM_WORLD,&A);
119: AssembleMatrix(PETSC_COMM_WORLD,&B);
120: MatShellSetOperation(B,MATOP_GET_SUBMATRIX,PETSC_NULL);
121: MatShellSetOperation(B,MATOP_GET_SUBMATRICES,PETSC_NULL);
122: MatGetOwnershipRange(A,&ms,PETSC_NULL);
124: idxrow[0] = ms+1;
125: idxrow[1] = ms+2;
126: idxrow[2] = ms+4;
127: ISCreateGeneral(PETSC_COMM_WORLD,3,idxrow,PETSC_USE_POINTER,&isrow);
129: idxcol[0] = ms+1;
130: idxcol[1] = ms+2;
131: idxcol[2] = ms+4;
132: idxcol[3] = ms+5;
133: ISCreateGeneral(PETSC_COMM_WORLD,4,idxcol,PETSC_USE_POINTER,&iscol);
135: MatGetSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,&Asub);
136: MatGetSubMatrix(B,isrow,iscol,MAT_INITIAL_MATRIX,&Bsub);
138: MatGetVecs(Asub,&right,&left);
139: VecDuplicate(right,&X);
140: VecDuplicate(right,&X1);
141: VecDuplicate(left,&Y);
142: VecDuplicate(left,&Y1);
144: PetscOptionsGetBool(PETSC_NULL,"-random",&random,PETSC_NULL);
145: if (random) {
146: VecSetRandom(right,PETSC_NULL);
147: VecSetRandom(left,PETSC_NULL);
148: VecSetRandom(X,PETSC_NULL);
149: VecSetRandom(Y,PETSC_NULL);
150: VecSetRandom(X1,PETSC_NULL);
151: VecSetRandom(Y1,PETSC_NULL);
152: } else {
153: VecSet(right,1.0);
154: VecSet(left,2.0);
155: VecSet(X,3.0);
156: VecSet(Y,4.0);
157: VecSet(X1,3.0);
158: VecSet(Y1,4.0);
159: }
160: CheckMatrices(Asub,Bsub,left,right,X,Y,X1,Y1);
161: ISDestroy(&isrow);
162: ISDestroy(&iscol);
163: MatDestroy(&A);
164: MatDestroy(&B);
165: MatDestroy(&Asub);
166: MatDestroy(&Bsub);
167: VecDestroy(&left);
168: VecDestroy(&right);
169: VecDestroy(&X);
170: VecDestroy(&Y);
171: VecDestroy(&X1);
172: VecDestroy(&Y1);
173: PetscFinalize();
174: return 0;
175: }