xref: /petsc/src/mat/tests/ex75.c (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce) 
1 
2 static char help[] = "Tests various routines in MatMPISBAIJ format.\n";
3 
4 #include <petscmat.h>
5 
6 int main(int argc,char **args)
7 {
8   Vec            x,y,u,s1,s2;
9   Mat            A,sA,sB;
10   PetscRandom    rctx;
11   PetscReal      r1,r2,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
12   PetscScalar    one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1;
13   PetscInt       n,col[3],n1,block,row,i,j,i2,j2,Ii,J,rstart,rend,bs=1,mbs=16,d_nz=3,o_nz=3,prob=1;
14   PetscErrorCode ierr;
15   PetscMPIInt    size,rank;
16   PetscBool      flg;
17 
18   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
19   ierr = PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);CHKERRQ(ierr);
20   ierr = PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);CHKERRQ(ierr);
21   ierr = PetscOptionsGetInt(NULL,NULL,"-prob",&prob,NULL);CHKERRQ(ierr);
22 
23   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
24   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
25 
26   /* Create a BAIJ matrix A */
27   n = mbs*bs;
28   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
29   ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
30   ierr = MatSetType(A,MATBAIJ);CHKERRQ(ierr);
31   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
32   ierr = MatMPIBAIJSetPreallocation(A,bs,d_nz,NULL,o_nz,NULL);CHKERRQ(ierr);
33   ierr = MatSeqBAIJSetPreallocation(A,bs,d_nz,NULL);CHKERRQ(ierr);
34   ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
35 
36   if (bs == 1) {
37     if (prob == 1) { /* tridiagonal matrix */
38       value[0] = -1.0; value[1] = 0.0; value[2] = -1.0;
39       for (i=1; i<n-1; i++) {
40         col[0] = i-1; col[1] = i; col[2] = i+1;
41         ierr   = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
42       }
43       i       = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
44       value[0]= 0.1; value[1]=-1.0; value[2]=0.0;
45       ierr    = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
46 
47       i        = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
48       value[0] = 0.0; value[1] = -1.0; value[2]=0.1;
49       ierr     = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
50     } else if (prob ==2) { /* matrix for the five point stencil */
51       n1 = (int) PetscSqrtReal((PetscReal)n);
52       for (i=0; i<n1; i++) {
53         for (j=0; j<n1; j++) {
54           Ii = j + n1*i;
55           if (i>0)    {J = Ii - n1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);}
56           if (i<n1-1) {J = Ii + n1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);}
57           if (j>0)    {J = Ii - 1;  ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);}
58           if (j<n1-1) {J = Ii + 1;  ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);}
59           ierr = MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);CHKERRQ(ierr);
60         }
61       }
62     }
63     /* end of if (bs == 1) */
64   } else {  /* bs > 1 */
65     for (block=0; block<n/bs; block++) {
66       /* diagonal blocks */
67       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
68       for (i=1+block*bs; i<bs-1+block*bs; i++) {
69         col[0] = i-1; col[1] = i; col[2] = i+1;
70         ierr   = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
71       }
72       i       = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
73       value[0]=-1.0; value[1]=4.0;
74       ierr    = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
75 
76       i       = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
77       value[0]=4.0; value[1] = -1.0;
78       ierr    = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
79     }
80     /* off-diagonal blocks */
81     value[0]=-1.0;
82     for (i=0; i<(n/bs-1)*bs; i++) {
83       col[0]=i+bs;
84       ierr  = MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);CHKERRQ(ierr);
85       col[0]=i; row=i+bs;
86       ierr  = MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr);
87     }
88   }
89   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
90   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
91   ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
92 
93   /* Get SBAIJ matrix sA from A */
94   ierr = MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr);
95 
96   /* Test MatGetSize(), MatGetLocalSize() */
97   ierr = MatGetSize(sA, &i,&j);CHKERRQ(ierr);
98   ierr = MatGetSize(A, &i2,&j2);CHKERRQ(ierr);
99   i   -= i2; j -= j2;
100   if (i || j) {
101     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);CHKERRQ(ierr);
102     ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
103   }
104 
105   ierr = MatGetLocalSize(sA, &i,&j);CHKERRQ(ierr);
106   ierr = MatGetLocalSize(A, &i2,&j2);CHKERRQ(ierr);
107   i2  -= i; j2 -= j;
108   if (i2 || j2) {
109     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);CHKERRQ(ierr);
110     ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
111   }
112 
113   /* vectors */
114   /*--------------------*/
115   /* i is obtained from MatGetLocalSize() */
116   ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
117   ierr = VecSetSizes(x,i,PETSC_DECIDE);CHKERRQ(ierr);
118   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
119   ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
120   ierr = VecDuplicate(x,&u);CHKERRQ(ierr);
121   ierr = VecDuplicate(x,&s1);CHKERRQ(ierr);
122   ierr = VecDuplicate(x,&s2);CHKERRQ(ierr);
123 
124   ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
125   ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
126   ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
127   ierr = VecSet(u,one);CHKERRQ(ierr);
128 
129   /* Test MatNorm() */
130   ierr  = MatNorm(A,NORM_FROBENIUS,&r1);CHKERRQ(ierr);
131   ierr  = MatNorm(sA,NORM_FROBENIUS,&r2);CHKERRQ(ierr);
132   rnorm = PetscAbsReal(r1-r2)/r2;
133   if (rnorm > tol && !rank) {
134     ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr);
135   }
136   ierr  = MatNorm(A,NORM_INFINITY,&r1);CHKERRQ(ierr);
137   ierr  = MatNorm(sA,NORM_INFINITY,&r2);CHKERRQ(ierr);
138   rnorm = PetscAbsReal(r1-r2)/r2;
139   if (rnorm > tol && !rank) {
140     ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr);
141   }
142   ierr  = MatNorm(A,NORM_1,&r1);CHKERRQ(ierr);
143   ierr  = MatNorm(sA,NORM_1,&r2);CHKERRQ(ierr);
144   rnorm = PetscAbsReal(r1-r2)/r2;
145   if (rnorm > tol && !rank) {
146     ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr);
147   }
148 
149   /* Test MatGetOwnershipRange() */
150   ierr = MatGetOwnershipRange(sA,&rstart,&rend);CHKERRQ(ierr);
151   ierr = MatGetOwnershipRange(A,&i2,&j2);CHKERRQ(ierr);
152   i2  -= rstart; j2 -= rend;
153   if (i2 || j2) {
154     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);CHKERRQ(ierr);
155     ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
156   }
157 
158   /* Test MatDiagonalScale() */
159   ierr = MatDiagonalScale(A,x,x);CHKERRQ(ierr);
160   ierr = MatDiagonalScale(sA,x,x);CHKERRQ(ierr);
161   ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr);
162   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
163 
164   /* Test MatGetDiagonal(), MatScale() */
165   ierr = MatGetDiagonal(A,s1);CHKERRQ(ierr);
166   ierr = MatGetDiagonal(sA,s2);CHKERRQ(ierr);
167   ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
168   ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
169   r1  -= r2;
170   if (r1<-tol || r1>tol) {
171     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);CHKERRQ(ierr);
172     ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
173   }
174 
175   ierr = MatScale(A,alpha);CHKERRQ(ierr);
176   ierr = MatScale(sA,alpha);CHKERRQ(ierr);
177 
178   /* Test MatGetRowMaxAbs() */
179   ierr = MatGetRowMaxAbs(A,s1,NULL);CHKERRQ(ierr);
180   ierr = MatGetRowMaxAbs(sA,s2,NULL);CHKERRQ(ierr);
181 
182   ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
183   ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
184   r1  -= r2;
185   if (r1<-tol || r1>tol) {
186     ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");CHKERRQ(ierr);
187   }
188 
189   /* Test MatMult(), MatMultAdd() */
190   ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr);
191   if (!flg) {
192     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);CHKERRQ(ierr);
193     ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
194   }
195 
196   ierr = MatMultAddEqual(A,sA,10,&flg);CHKERRQ(ierr);
197   if (!flg) {
198     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);CHKERRQ(ierr);
199     ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
200   }
201 
202   /* Test MatMultTranspose(), MatMultTransposeAdd() */
203   for (i=0; i<10; i++) {
204     ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
205     ierr = MatMultTranspose(A,x,s1);CHKERRQ(ierr);
206     ierr = MatMultTranspose(sA,x,s2);CHKERRQ(ierr);
207     ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
208     ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
209     r1  -= r2;
210     if (r1<-tol || r1>tol) {
211       ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);CHKERRQ(ierr);
212       ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
213     }
214   }
215   for (i=0; i<10; i++) {
216     ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
217     ierr = VecSetRandom(y,rctx);CHKERRQ(ierr);
218     ierr = MatMultTransposeAdd(A,x,y,s1);CHKERRQ(ierr);
219     ierr = MatMultTransposeAdd(sA,x,y,s2);CHKERRQ(ierr);
220     ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
221     ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
222     r1  -= r2;
223     if (r1<-tol || r1>tol) {
224       ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);CHKERRQ(ierr);
225       ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
226     }
227   }
228 
229   /* Test MatDuplicate() */
230   ierr = MatDuplicate(sA,MAT_COPY_VALUES,&sB);CHKERRQ(ierr);
231   ierr = MatEqual(sA,sB,&flg);CHKERRQ(ierr);
232   if (!flg) {
233     ierr = PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");CHKERRQ(ierr);
234   }
235   ierr = MatMultEqual(sA,sB,5,&flg);CHKERRQ(ierr);
236   if (!flg) {
237     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);CHKERRQ(ierr);
238     ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
239   }
240   ierr = MatMultAddEqual(sA,sB,5,&flg);CHKERRQ(ierr);
241   if (!flg) {
242     ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);CHKERRQ(ierr);
243     ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
244   }
245   ierr = MatDestroy(&sB);CHKERRQ(ierr);
246   ierr = VecDestroy(&u);CHKERRQ(ierr);
247   ierr = VecDestroy(&x);CHKERRQ(ierr);
248   ierr = VecDestroy(&y);CHKERRQ(ierr);
249   ierr = VecDestroy(&s1);CHKERRQ(ierr);
250   ierr = VecDestroy(&s2);CHKERRQ(ierr);
251   ierr = MatDestroy(&sA);CHKERRQ(ierr);
252   ierr = MatDestroy(&A);CHKERRQ(ierr);
253   ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
254   ierr = PetscFinalize();
255   return ierr;
256 }
257 
258 /*TEST
259 
260    test:
261       nsize: {{1 3}}
262       args: -bs {{1 2 3  5  7 8}} -mat_ignore_lower_triangular -prob {{1 2}}
263 
264 TEST*/
265