xref: /petsc/src/mat/impls/baij/seq/baij2.c (revision 83e2fdc756c13f964f3ed5aaff42e3626c5b72bc) 
1 #ifdef PETSC_RCS_HEADER
2 static char vcid[] = "$Id: baij2.c,v 1.15 1997/07/09 20:55:07 balay Exp bsmith $";
3 #endif
4 
5 #include "src/mat/impls/baij/seq/baij.h"
6 #include "petsc.h"
7 #include "src/inline/bitarray.h"
8 
9 #undef __FUNC__
10 #define __FUNC__ "MatIncreaseOverlap_SeqBAIJ" /* ADIC Ignore */
11 int MatIncreaseOverlap_SeqBAIJ(Mat A,int is_max,IS *is,int ov)
12 {
13   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data;
14   int         row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val, ival;
15   int         start, end, *ai, *aj,bs,*nidx2;
16   char        *table;
17 
18   m     = a->mbs;
19   ai    = a->i;
20   aj    = a->j;
21   bs    = a->bs;
22 
23   if (ov < 0)  SETERRQ(1,0,"Negative overlap specified");
24 
25   table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table);
26   nidx  = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx);
27   nidx2 = (int *)PetscMalloc((a->m+1)*sizeof(int)); CHKPTRQ(nidx2);
28 
29   for ( i=0; i<is_max; i++ ) {
30     /* Initialise the two local arrays */
31     isz  = 0;
32     PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char));
33 
34     /* Extract the indices, assume there can be duplicate entries */
35     ierr = ISGetIndices(is[i],&idx);  CHKERRQ(ierr);
36     ierr = ISGetSize(is[i],&n);  CHKERRQ(ierr);
37 
38     /* Enter these into the temp arrays i.e mark table[row], enter row into new index */
39     for ( j=0; j<n ; ++j){
40       ival = idx[j]/bs; /* convert the indices into block indices */
41       if (ival>m) SETERRQ(1,0,"index greater than mat-dim");
42       if(!BT_LOOKUP(table, ival)) { nidx[isz++] = ival;}
43     }
44     ierr = ISRestoreIndices(is[i],&idx);  CHKERRQ(ierr);
45     ierr = ISDestroy(is[i]); CHKERRQ(ierr);
46 
47     k = 0;
48     for ( j=0; j<ov; j++){ /* for each overlap*/
49       n = isz;
50       for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */
51         row   = nidx[k];
52         start = ai[row];
53         end   = ai[row+1];
54         for ( l = start; l<end ; l++){
55           val = aj[l];
56           if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;}
57         }
58       }
59     }
60     /* expand the Index Set */
61     for (j=0; j<isz; j++ ) {
62       for (k=0; k<bs; k++ )
63         nidx2[j*bs+k] = nidx[j]*bs+k;
64     }
65     ierr = ISCreateGeneral(PETSC_COMM_SELF, isz*bs, nidx2, (is+i)); CHKERRQ(ierr);
66   }
67   PetscFree(table);
68   PetscFree(nidx);
69   PetscFree(nidx2);
70   return 0;
71 }
72 
73 #undef __FUNC__
74 #define __FUNC__ "MatGetSubMatrix_SeqBAIJ_Private"
75 int MatGetSubMatrix_SeqBAIJ_Private(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B)
76 {
77   Mat_SeqBAIJ  *a = (Mat_SeqBAIJ *) A->data,*c;
78   int          nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->nbs,*lens;
79   int          row,mat_i,*mat_j,tcol,*mat_ilen;
80   int          *irow, *icol, nrows, ncols,*ssmap,bs=a->bs, bs2=a->bs2;
81   int          *aj = a->j, *ai = a->i;
82   Scalar       *mat_a;
83   Mat          C;
84 
85   ierr = ISSorted(iscol,(PetscTruth*)&i);
86   if (!i) SETERRQ(1,0,"IS is not sorted");
87 
88   ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr);
89   ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr);
90   ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr);
91   ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr);
92 
93   smap  = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap);
94   ssmap = smap;
95   lens  = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens);
96   PetscMemzero(smap,oldcols*sizeof(int));
97   for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1;
98   /* determine lens of each row */
99   for (i=0; i<nrows; i++) {
100     kstart  = ai[irow[i]];
101     kend    = kstart + a->ilen[irow[i]];
102     lens[i] = 0;
103       for ( k=kstart; k<kend; k++ ) {
104         if (ssmap[aj[k]]) {
105           lens[i]++;
106         }
107       }
108     }
109   /* Create and fill new matrix */
110   if (scall == MAT_REUSE_MATRIX) {
111     c = (Mat_SeqBAIJ *)((*B)->data);
112 
113     if (c->mbs!=nrows || c->nbs!=ncols || c->bs!=bs)
114       SETERRQ(1,0,"");
115     if (PetscMemcmp(c->ilen,lens, c->mbs *sizeof(int))) {
116       SETERRQ(1,0,"Cannot reuse matrix. wrong no of nonzeros");
117     }
118     PetscMemzero(c->ilen,c->mbs*sizeof(int));
119     C = *B;
120     }
121   else {
122     ierr = MatCreateSeqBAIJ(A->comm,bs,nrows*bs,ncols*bs,0,lens,&C);CHKERRQ(ierr);
123   }
124   c = (Mat_SeqBAIJ *)(C->data);
125   for (i=0; i<nrows; i++) {
126     row    = irow[i];
127     nznew  = 0;
128     kstart = ai[row];
129     kend   = kstart + a->ilen[row];
130     mat_i  = c->i[i];
131     mat_j  = c->j + mat_i;
132     mat_a  = c->a + mat_i*bs2;
133     mat_ilen = c->ilen + i;
134     for ( k=kstart; k<kend; k++ ) {
135       if ((tcol=ssmap[a->j[k]])) {
136         *mat_j++ = tcol - 1;
137         PetscMemcpy(mat_a,a->a+k*bs2,bs2*sizeof(Scalar)); mat_a+=bs2;
138         (*mat_ilen)++;
139 
140       }
141     }
142   }
143 
144   /* Free work space */
145   ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr);
146   PetscFree(smap); PetscFree(lens);
147   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
148   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
149 
150   ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr);
151   *B = C;
152   return 0;
153 }
154 
155 #undef __FUNC__
156 #define __FUNC__ "MatGetSubMatrix_SeqBAIJ"
157 int MatGetSubMatrix_SeqBAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B)
158 {
159   Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data;
160   IS          is1,is2;
161   int         *vary,*iary,*irow,*icol,nrows,ncols,i,ierr,bs=a->bs,count;
162 
163   ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr);
164   ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr);
165   ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr);
166   ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr);
167 
168   /* Verify if the indices corespond to each elementin a block
169    and form the IS with compressed IS */
170   vary = (int *) PetscMalloc(2*(a->mbs+1)*sizeof(int)); CHKPTRQ(vary);
171   iary = vary + a->mbs;
172   PetscMemzero(vary,(a->mbs)*sizeof(int));
173   for ( i=0; i<nrows; i++) vary[irow[i]/bs]++;
174   count = 0;
175   for (i=0; i<a->mbs; i++) {
176     if (vary[i]!=0 && vary[i]!=bs) SETERRA(1,0,"MatGetSubmatrices_SeqBAIJ:");
177     if (vary[i]==bs) iary[count++] = i;
178   }
179   ierr = ISCreateGeneral(PETSC_COMM_SELF, count, iary,&is1); CHKERRQ(ierr);
180 
181   PetscMemzero(vary,(a->mbs)*sizeof(int));
182   for ( i=0; i<ncols; i++) vary[icol[i]/bs]++;
183   count = 0;
184   for (i=0; i<a->mbs; i++) {
185     if (vary[i]!=0 && vary[i]!=bs) SETERRA(1,0,"MatGetSubmatrices_SeqBAIJ:");
186     if (vary[i]==bs) iary[count++] = i;
187   }
188   ierr = ISCreateGeneral(PETSC_COMM_SELF, count, iary,&is2); CHKERRQ(ierr);
189   ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr);
190   ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr);
191   PetscFree(vary);
192 
193   ierr = MatGetSubMatrix_SeqBAIJ_Private(A,is1,is2,scall,B); CHKERRQ(ierr);
194   ISDestroy(is1);
195   ISDestroy(is2);
196   return 0;
197 }
198 
199 extern int MatGetSubMatrix_SeqBAIJ(Mat,IS,IS,MatGetSubMatrixCall,Mat*);
200 
201 #undef __FUNC__
202 #define __FUNC__ "MatGetSubMatrices_SeqBAIJ"
203 int MatGetSubMatrices_SeqBAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,
204                                     Mat **B)
205 {
206   int ierr,i;
207 
208   if (scall == MAT_INITIAL_MATRIX) {
209     *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B);
210   }
211 
212   for ( i=0; i<n; i++ ) {
213     ierr = MatGetSubMatrix_SeqBAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr);
214   }
215   return 0;
216 }
217 
218 
219 
220 
221 
222 
223 
224