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