xref: /petsc/src/mat/impls/sbaij/seq/aijsbaij.c (revision 9c30b7d2697335155d7490a7e085415ee7b4a02a) 
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_SeqSBAI_SeqAIJ"
9 int MatConvert_SeqSBAIJ_SeqAIJ(Mat A,const MatType newtype,Mat *newmat)
10 {
11   Mat          B;
12   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
13   Mat_SeqAIJ   *b;
14   int          ierr,*ai=a->i,*aj=a->j,m=A->m,n=A->n,i,j,k,*bi,*bj,
15                *rowlengths,nz,*rowstart,itmp;
16   int          bs=a->bs,bs2=bs*bs,mbs=A->m/bs;
17   PetscScalar  *av,*bv;
18 
19   PetscFunctionBegin;
20 
21   /* compute rowlengths of newmat */
22   ierr = PetscMalloc((2*m+1)*sizeof(int),&rowlengths);CHKERRQ(ierr);
23   rowstart = rowlengths + m;
24 
25   for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
26   aj = a->j;
27   k = 0;
28   for (i=0; i<mbs; i++) {
29     nz = ai[i+1] - ai[i];
30     aj++; /* skip diagonal */
31     for (j=1; j<nz; j++) { /* no. of lower triangular blocks */
32       rowlengths[(*aj)*bs]++; aj++;
33     }
34     rowlengths[k] += nz;   /* no. of upper triangular blocks */
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(A->comm,m,n,m,n,&B);CHKERRQ(ierr);
44   ierr = MatSetType(B,newtype);CHKERRQ(ierr);
45   ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr);
46   ierr = MatSetOption(B,MAT_COLUMN_ORIENTED);CHKERRQ(ierr);
47   ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr);
48   ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr);
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 - mbs) SETERRQ2(1,"bi[mbs]: %d != 2*a->nz-mbs: %d\n",bi[mbs],2*a->nz - mbs);
65 
66   /* set b->j and b->a */
67   aj = a->j; av = a->a;
68   for (i=0; i<mbs; i++) {
69     /* diagonal block */
70     for (j=0; j<bs; j++){   /* row i*bs+j */
71       itmp = i*bs+j;
72       for (k=0; k<bs; k++){ /* col i*bs+k */
73         *(bj + rowstart[itmp]) = (*aj)*bs+k;
74         *(bv + rowstart[itmp]) = *(av+k*bs+j);
75         rowstart[itmp]++;
76       }
77     }
78     aj++; av += bs2;
79 
80     nz = ai[i+1] - ai[i] -1;
81     while (nz--){
82       /* lower triangular blocks */
83       for (j=0; j<bs; j++){   /* row (*aj)*bs+j */
84         itmp = (*aj)*bs+j;
85         for (k=0; k<bs; k++){ /* col i*bs+k */
86           *(bj + rowstart[itmp]) = i*bs+k;
87           *(bv + rowstart[itmp]) = *(av+k*bs+j);
88           rowstart[itmp]++;
89         }
90       }
91       /* upper triangular blocks */
92       for (j=0; j<bs; j++){   /* row i*bs+j */
93         itmp = i*bs+j;
94         for (k=0; k<bs; k++){ /* col (*aj)*bs+k */
95           *(bj + rowstart[itmp]) = (*aj)*bs+k;
96           *(bv + rowstart[itmp]) = *(av+k*bs+j);
97           rowstart[itmp]++;
98         }
99       }
100       aj++; av += bs2;
101     }
102   }
103   ierr = PetscFree(rowlengths);CHKERRQ(ierr);
104   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
105   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
106 
107   /* Fake support for "inplace" convert. */
108   if (*newmat == A) {
109     ierr = MatDestroy(A);CHKERRQ(ierr);
110   }
111   *newmat = B;
112   PetscFunctionReturn(0);
113 }
114 #undef __FUNCT__
115 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ"
116 int MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,Mat *newmat) {
117   Mat          B;
118   Mat_SeqAIJ   *a = (Mat_SeqAIJ*)A->data;
119   Mat_SeqSBAIJ *b;
120   int          ierr,*ai=a->i,*aj,m=A->M,n=A->N,i,j,
121                *bi,*bj,*rowlengths;
122   PetscScalar  *av,*bv;
123 
124   PetscFunctionBegin;
125   if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
126   ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */
127 
128   ierr = PetscMalloc(m*sizeof(int),&rowlengths);CHKERRQ(ierr);
129   for (i=0; i<m; i++) {
130     rowlengths[i] = ai[i+1] - a->diag[i];
131   }
132   ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr);
133   ierr = MatSetType(B,newtype);CHKERRQ(ierr);
134   ierr = MatSeqSBAIJSetPreallocation(B,1,0,rowlengths);CHKERRQ(ierr);
135 
136   ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr);
137   ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr);
138   ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr);
139 
140   b  = (Mat_SeqSBAIJ*)(B->data);
141   bi = b->i;
142   bj = b->j;
143   bv = b->a;
144 
145   bi[0] = 0;
146   for (i=0; i<m; i++) {
147     aj = a->j + a->diag[i];
148     av = a->a + a->diag[i];
149     for (j=0; j<rowlengths[i]; j++){
150       *bj = *aj; bj++; aj++;
151       *bv = *av; bv++; av++;
152     }
153     bi[i+1]    = bi[i] + rowlengths[i];
154     b->ilen[i] = rowlengths[i];
155   }
156 
157   ierr = PetscFree(rowlengths);CHKERRQ(ierr);
158   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
159   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
160 
161   /* Fake support for "inplace" convert. */
162   if (*newmat == A) {
163     ierr = MatDestroy(A);CHKERRQ(ierr);
164   }
165   *newmat = B;
166 
167   PetscFunctionReturn(0);
168 }
169 EXTERN_C_END
170 
171 EXTERN_C_BEGIN
172 #undef __FUNCT__
173 #define __FUNCT__ "MatConvert_SeqSBAI_SeqBAIJ"
174 int MatConvert_SeqSBAIJ_SeqBAIJ(Mat A,const MatType newtype,Mat *newmat)
175 {
176   Mat          B;
177   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
178   Mat_SeqBAIJ  *b;
179   int          ierr,*ai=a->i,*aj=a->j,m=A->m,n=A->n,i,k,*bi,*bj,
180                *browlengths,nz,*browstart,itmp;
181   int          bs=a->bs,bs2=bs*bs,mbs=m/bs;
182   PetscScalar  *av,*bv;
183 
184   PetscFunctionBegin;
185   /* compute browlengths of newmat */
186   ierr = PetscMalloc(2*mbs*sizeof(int),&browlengths);CHKERRQ(ierr);
187   browstart = browlengths + mbs;
188   for (i=0; i<mbs; i++) browlengths[i] = 0;
189   aj = a->j;
190   for (i=0; i<mbs; i++) {
191     nz = ai[i+1] - ai[i];
192     aj++; /* skip diagonal */
193     for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
194       browlengths[*aj]++; aj++;
195     }
196     browlengths[i] += nz;   /* no. of upper triangular blocks */
197   }
198 
199   ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr);
200   ierr = MatSetType(B,newtype);CHKERRQ(ierr);
201   ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr);
202   ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr);
203   ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr);
204   ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr);
205 
206   b  = (Mat_SeqBAIJ*)(B->data);
207   bi = b->i;
208   bj = b->j;
209   bv = b->a;
210 
211   /* set b->i */
212   bi[0] = 0;
213   for (i=0; i<mbs; i++){
214     b->ilen[i]    = browlengths[i];
215     bi[i+1]       = bi[i] + browlengths[i];
216     browstart[i]  = bi[i];
217   }
218   if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(1,"bi[mbs]: %d != 2*a->nz - mbs: %d\n",bi[mbs],2*a->nz - mbs);
219 
220   /* set b->j and b->a */
221   aj = a->j; av = a->a;
222   for (i=0; i<mbs; i++) {
223     /* diagonal block */
224     *(bj + browstart[i]) = *aj; aj++;
225     itmp = bs2*browstart[i];
226     for (k=0; k<bs2; k++){
227       *(bv + itmp + k) = *av; av++;
228     }
229     browstart[i]++;
230 
231     nz = ai[i+1] - ai[i] -1;
232     while (nz--){
233       /* lower triangular blocks */
234       *(bj + browstart[*aj]) = i;
235       itmp = bs2*browstart[*aj];
236       for (k=0; k<bs2; k++){
237         *(bv + itmp + k) = *(av + k);
238       }
239       browstart[*aj]++;
240 
241       /* upper triangular blocks */
242       *(bj + browstart[i]) = *aj; aj++;
243       itmp = bs2*browstart[i];
244       for (k=0; k<bs2; k++){
245         *(bv + itmp + k) = *av; av++;
246       }
247       browstart[i]++;
248     }
249   }
250   ierr = PetscFree(browlengths);CHKERRQ(ierr);
251   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
252   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
253 
254   /* Fake support for "inplace" convert. */
255   if (*newmat == A) {
256     ierr = MatDestroy(A);CHKERRQ(ierr);
257   }
258   *newmat = B;
259   PetscFunctionReturn(0);
260 }
261 #undef __FUNCT__
262 #define __FUNCT__ "MatConvert_SeqBAIJ_SeqSBAIJ"
263 int MatConvert_SeqBAIJ_SeqSBAIJ(Mat A,const MatType newtype,Mat *newmat)
264 {
265   Mat          B;
266   Mat_SeqBAIJ  *a = (Mat_SeqBAIJ*)A->data;
267   Mat_SeqSBAIJ *b;
268   int          ierr,*ai=a->i,*aj,m=A->m,n=A->n,i,j,k,
269                *bi,*bj,*browlengths;
270   int          bs=a->bs,bs2=bs*bs,mbs=m/bs;
271   PetscScalar  *av,*bv;
272 
273   PetscFunctionBegin;
274   if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
275   ierr = MatMissingDiagonal_SeqBAIJ(A);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */
276 
277   ierr = PetscMalloc(mbs*sizeof(int),&browlengths);CHKERRQ(ierr);
278   for (i=0; i<mbs; i++) {
279     browlengths[i] = ai[i+1] - a->diag[i];
280   }
281 
282   ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr);
283   ierr = MatSetType(B,newtype);CHKERRQ(ierr);
284   ierr = MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr);
285   ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr);
286   ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr);
287   ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr);
288 
289   b  = (Mat_SeqSBAIJ*)(B->data);
290   bi = b->i;
291   bj = b->j;
292   bv = b->a;
293 
294   bi[0] = 0;
295   for (i=0; i<mbs; i++) {
296     aj = a->j + a->diag[i];
297     av = a->a + (a->diag[i])*bs2;
298     for (j=0; j<browlengths[i]; j++){
299       *bj = *aj; bj++; aj++;
300       for (k=0; k<bs2; k++){
301         *bv = *av; bv++; av++;
302       }
303     }
304     bi[i+1]    = bi[i] + browlengths[i];
305     b->ilen[i] = browlengths[i];
306   }
307   ierr = PetscFree(browlengths);CHKERRQ(ierr);
308   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
309   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
310 
311   /* Fake support for "inplace" convert. */
312   if (*newmat == A) {
313     ierr = MatDestroy(A);CHKERRQ(ierr);
314   }
315   *newmat = B;
316 
317   PetscFunctionReturn(0);
318 }
319 EXTERN_C_END
320