xref: /petsc/src/mat/impls/sbaij/seq/aijsbaij.c (revision c5b2ea42cd4070151795b0df3c8c47a133418b2b) 
1 
2 #include <../src/mat/impls/aij/seq/aij.h>
3 #include <../src/mat/impls/baij/seq/baij.h>
4 #include <../src/mat/impls/sbaij/seq/sbaij.h>
5 
6 EXTERN_C_BEGIN
7 #undef __FUNCT__
8 #define __FUNCT__ "MatConvert_SeqSBAIJ_SeqAIJ"
9 PetscErrorCode  MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
10 {
11   Mat            B;
12   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data;
13   Mat_SeqAIJ     *b;
14   PetscErrorCode ierr;
15   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp;
16   PetscInt       bs=A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs,diagcnt=0;
17   MatScalar      *av,*bv;
18 
19   PetscFunctionBegin;
20   /* compute rowlengths of newmat */
21   ierr = PetscMalloc2(m,PetscInt,&rowlengths,m+1,PetscInt,&rowstart);CHKERRQ(ierr);
22 
23   for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
24   aj = a->j;
25   k = 0;
26   for (i=0; i<mbs; i++) {
27     nz = ai[i+1] - ai[i];
28     if (nz) {
29       rowlengths[k] += nz;   /* no. of upper triangular blocks */
30       if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */
31       for (j=0; j<nz; j++) { /* no. of lower triangular blocks */
32 	rowlengths[(*aj)*bs]++; aj++;
33       }
34     }
35     rowlengths[k] *= bs;
36     for (j=1; j<bs; j++) {
37       rowlengths[k+j] = rowlengths[k];
38     }
39     k += bs;
40     /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */
41   }
42 
43   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
44   ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr);
45   ierr = MatSetType(B,newtype);CHKERRQ(ierr);
46   ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr);
47   ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr);
48   B->rmap->bs = A->rmap->bs;
49 
50   b  = (Mat_SeqAIJ*)(B->data);
51   bi = b->i;
52   bj = b->j;
53   bv = b->a;
54 
55   /* set b->i */
56   bi[0] = 0; rowstart[0] = 0;
57   for (i=0; i<mbs; i++){
58     for (j=0; j<bs; j++){
59       b->ilen[i*bs+j]    = rowlengths[i*bs];
60       rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs];
61     }
62     bi[i+1]     = bi[i] + rowlengths[i*bs]/bs;
63   }
64   if (bi[mbs] != 2*a->nz - diagcnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-diagcnt: %D\n",bi[mbs],2*a->nz - diagcnt);
65 
66   /* set b->j and b->a */
67   aj = a->j; av = a->a;
68   for (i=0; i<mbs; i++) {
69     nz = ai[i+1] - ai[i];
70     /* diagonal block */
71     if (nz && *aj == i) {
72       nz--;
73       for (j=0; j<bs; j++){   /* row i*bs+j */
74 	itmp = i*bs+j;
75 	for (k=0; k<bs; k++){ /* col i*bs+k */
76 	  *(bj + rowstart[itmp]) = (*aj)*bs+k;
77 	  *(bv + rowstart[itmp]) = *(av+k*bs+j);
78 	  rowstart[itmp]++;
79 	}
80       }
81       aj++; av += bs2;
82     }
83 
84     while (nz--){
85       /* lower triangular blocks */
86       for (j=0; j<bs; j++){   /* row (*aj)*bs+j */
87         itmp = (*aj)*bs+j;
88         for (k=0; k<bs; k++){ /* col i*bs+k */
89           *(bj + rowstart[itmp]) = i*bs+k;
90           *(bv + rowstart[itmp]) = *(av+j*bs+k);
91           rowstart[itmp]++;
92         }
93       }
94       /* upper triangular blocks */
95       for (j=0; j<bs; j++){   /* row i*bs+j */
96         itmp = i*bs+j;
97         for (k=0; k<bs; k++){ /* col (*aj)*bs+k */
98           *(bj + rowstart[itmp]) = (*aj)*bs+k;
99           *(bv + rowstart[itmp]) = *(av+k*bs+j);
100           rowstart[itmp]++;
101         }
102       }
103       aj++; av += bs2;
104     }
105   }
106   ierr = PetscFree2(rowlengths,rowstart);CHKERRQ(ierr);
107   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
108   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
109 
110   if (reuse == MAT_REUSE_MATRIX) {
111     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
112   } else {
113     *newmat = B;
114   }
115   PetscFunctionReturn(0);
116 }
117 EXTERN_C_END
118 
119 EXTERN_C_BEGIN
120 #undef __FUNCT__
121 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ"
122 PetscErrorCode  MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat) {
123   Mat            B;
124   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
125   Mat_SeqSBAIJ   *b;
126   PetscErrorCode ierr;
127   PetscInt       *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths;
128   MatScalar      *av,*bv;
129 
130   PetscFunctionBegin;
131   if (!A->symmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
132   if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
133 
134   ierr = PetscMalloc(m*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr);
135   for (i=0; i<m; i++) {
136     rowlengths[i] = ai[i+1] - a->diag[i];
137   }
138   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
139   ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr);
140   ierr = MatSetType(B,newtype);CHKERRQ(ierr);
141   ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,1,0,rowlengths);CHKERRQ(ierr);
142 
143   ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr);
144 
145   b  = (Mat_SeqSBAIJ*)(B->data);
146   bi = b->i;
147   bj = b->j;
148   bv = b->a;
149 
150   bi[0] = 0;
151   for (i=0; i<m; i++) {
152     aj = a->j + a->diag[i];
153     av = a->a + a->diag[i];
154     for (j=0; j<rowlengths[i]; j++){
155       *bj = *aj; bj++; aj++;
156       *bv = *av; bv++; av++;
157     }
158     bi[i+1]    = bi[i] + rowlengths[i];
159     b->ilen[i] = rowlengths[i];
160   }
161 
162   ierr = PetscFree(rowlengths);CHKERRQ(ierr);
163   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
164   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
165   if (A->hermitian){
166     ierr = MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr);
167   }
168 
169   if (reuse == MAT_REUSE_MATRIX) {
170     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
171   } else {
172     *newmat = B;
173   }
174   PetscFunctionReturn(0);
175 }
176 EXTERN_C_END
177 
178 EXTERN_C_BEGIN
179 #undef __FUNCT__
180 #define __FUNCT__ "MatConvert_SeqSBAIJ_SeqBAIJ"
181 PetscErrorCode  MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
182 {
183   Mat            B;
184   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data;
185   Mat_SeqBAIJ    *b;
186   PetscErrorCode ierr;
187   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
188   PetscInt       bs=A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
189   MatScalar      *av,*bv;
190 
191   PetscFunctionBegin;
192   /* compute browlengths of newmat */
193   ierr = PetscMalloc2(mbs,PetscInt,&browlengths,mbs,PetscInt,&browstart);CHKERRQ(ierr);
194   for (i=0; i<mbs; i++) browlengths[i] = 0;
195   aj = a->j;
196   for (i=0; i<mbs; i++) {
197     nz = ai[i+1] - ai[i];
198     aj++; /* skip diagonal */
199     for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
200       browlengths[*aj]++; aj++;
201     }
202     browlengths[i] += nz;   /* no. of upper triangular blocks */
203   }
204 
205   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
206   ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr);
207   ierr = MatSetType(B,newtype);CHKERRQ(ierr);
208   ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr);
209   ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr);
210 
211   b  = (Mat_SeqBAIJ*)(B->data);
212   bi = b->i;
213   bj = b->j;
214   bv = b->a;
215 
216   /* set b->i */
217   bi[0] = 0;
218   for (i=0; i<mbs; i++){
219     b->ilen[i]   = browlengths[i];
220     bi[i+1]      = bi[i] + browlengths[i];
221     browstart[i] = bi[i];
222   }
223   if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz - mbs: %D\n",bi[mbs],2*a->nz - mbs);
224 
225   /* set b->j and b->a */
226   aj = a->j; av = a->a;
227   for (i=0; i<mbs; i++) {
228     /* diagonal block */
229     *(bj + browstart[i]) = *aj; aj++;
230     itmp = bs2*browstart[i];
231     for (k=0; k<bs2; k++){
232       *(bv + itmp + k) = *av; av++;
233     }
234     browstart[i]++;
235 
236     nz = ai[i+1] - ai[i] -1;
237     while (nz--){
238       /* lower triangular blocks - transpose blocks of A */
239       *(bj + browstart[*aj]) = i; /* block col index */
240       itmp = bs2*browstart[*aj];  /* row index */
241       for (col=0; col<bs; col++){
242         k = col;
243         for (row=0; row<bs; row++){
244           bv[itmp + col*bs+row] = av[k]; k+=bs;
245         }
246       }
247       browstart[*aj]++;
248 
249       /* upper triangular blocks */
250       *(bj + browstart[i]) = *aj; aj++;
251       itmp = bs2*browstart[i];
252       for (k=0; k<bs2; k++){
253         bv[itmp + k] = av[k];
254       }
255       av += bs2;
256       browstart[i]++;
257     }
258   }
259   ierr = PetscFree2(browlengths,browstart);CHKERRQ(ierr);
260   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
261   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
262 
263   if (reuse == MAT_REUSE_MATRIX) {
264     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
265   } else {
266     *newmat = B;
267   }
268   PetscFunctionReturn(0);
269 }
270 EXTERN_C_END
271 
272 EXTERN_C_BEGIN
273 #undef __FUNCT__
274 #define __FUNCT__ "MatConvert_SeqBAIJ_SeqSBAIJ"
275 PetscErrorCode  MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
276 {
277   Mat            B;
278   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
279   Mat_SeqSBAIJ   *b;
280   PetscErrorCode ierr;
281   PetscInt       *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
282   PetscInt       bs=A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
283   MatScalar      *av,*bv;
284   PetscBool      flg;
285 
286   PetscFunctionBegin;
287   if (!A->symmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
288   if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
289   ierr = MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */
290   if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd);
291 
292   ierr = PetscMalloc(mbs*sizeof(PetscInt),&browlengths);CHKERRQ(ierr);
293   for (i=0; i<mbs; i++) {
294     browlengths[i] = ai[i+1] - a->diag[i];
295   }
296 
297   ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
298   ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr);
299   ierr = MatSetType(B,newtype);CHKERRQ(ierr);
300   ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,browlengths);CHKERRQ(ierr);
301   ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr);
302 
303   b  = (Mat_SeqSBAIJ*)(B->data);
304   bi = b->i;
305   bj = b->j;
306   bv = b->a;
307 
308   bi[0] = 0;
309   for (i=0; i<mbs; i++) {
310     aj = a->j + a->diag[i];
311     av = a->a + (a->diag[i])*bs2;
312     for (j=0; j<browlengths[i]; j++){
313       *bj = *aj; bj++; aj++;
314       for (k=0; k<bs2; k++){
315         *bv = *av; bv++; av++;
316       }
317     }
318     bi[i+1]    = bi[i] + browlengths[i];
319     b->ilen[i] = browlengths[i];
320   }
321   ierr = PetscFree(browlengths);CHKERRQ(ierr);
322   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
323   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
324 
325   if (reuse == MAT_REUSE_MATRIX) {
326     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
327   } else {
328     *newmat = B;
329   }
330 
331   PetscFunctionReturn(0);
332 }
333 EXTERN_C_END
334