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