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