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