xref: /petsc/src/mat/impls/baij/seq/baij2.c (revision 3a40ed3dce77c081171d005ae1a6ff4bb9d13b6f) 
1 #ifdef PETSC_RCS_HEADER
2 static char vcid[] = "$Id: baij2.c,v 1.18 1997/10/01 21:30:01 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"
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   BT          table;
17 
18   PetscFunctionBegin;
19   m     = a->mbs;
20   ai    = a->i;
21   aj    = a->j;
22   bs    = a->bs;
23 
24   if (ov < 0)  SETERRQ(1,0,"Negative overlap specified");
25 
26   ierr  = BTCreate(m,table); CHKERRQ(ierr);
27   nidx  = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx);
28   nidx2 = (int *)PetscMalloc((a->m+1)*sizeof(int)); CHKPTRQ(nidx2);
29 
30   for ( i=0; i<is_max; i++ ) {
31     /* Initialise the two local arrays */
32     isz  = 0;
33     BTMemzero(m,table);
34 
35     /* Extract the indices, assume there can be duplicate entries */
36     ierr = ISGetIndices(is[i],&idx);  CHKERRQ(ierr);
37     ierr = ISGetSize(is[i],&n);  CHKERRQ(ierr);
38 
39     /* Enter these into the temp arrays i.e mark table[row], enter row into new index */
40     for ( j=0; j<n ; ++j){
41       ival = idx[j]/bs; /* convert the indices into block indices */
42       if (ival>m) SETERRQ(1,0,"index greater than mat-dim");
43       if(!BTLookupSet(table, ival)) { nidx[isz++] = ival;}
44     }
45     ierr = ISRestoreIndices(is[i],&idx);  CHKERRQ(ierr);
46     ierr = ISDestroy(is[i]); CHKERRQ(ierr);
47 
48     k = 0;
49     for ( j=0; j<ov; j++){ /* for each overlap*/
50       n = isz;
51       for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */
52         row   = nidx[k];
53         start = ai[row];
54         end   = ai[row+1];
55         for ( l = start; l<end ; l++){
56           val = aj[l];
57           if (!BTLookupSet(table,val)) {nidx[isz++] = val;}
58         }
59       }
60     }
61     /* expand the Index Set */
62     for (j=0; j<isz; j++ ) {
63       for (k=0; k<bs; k++ )
64         nidx2[j*bs+k] = nidx[j]*bs+k;
65     }
66     ierr = ISCreateGeneral(PETSC_COMM_SELF, isz*bs, nidx2, (is+i)); CHKERRQ(ierr);
67   }
68   BTDestroy(table);
69   PetscFree(nidx);
70   PetscFree(nidx2);
71   PetscFunctionReturn(0);
72 }
73 
74 #undef __FUNC__
75 #define __FUNC__ "MatGetSubMatrix_SeqBAIJ_Private"
76 int MatGetSubMatrix_SeqBAIJ_Private(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B)
77 {
78   Mat_SeqBAIJ  *a = (Mat_SeqBAIJ *) A->data,*c;
79   int          nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->nbs,*lens;
80   int          row,mat_i,*mat_j,tcol,*mat_ilen;
81   int          *irow, *icol, nrows, ncols,*ssmap,bs=a->bs, bs2=a->bs2;
82   int          *aj = a->j, *ai = a->i;
83   Scalar       *mat_a;
84   Mat          C;
85 
86   PetscFunctionBegin;
87   ierr = ISSorted(iscol,(PetscTruth*)&i);
88   if (!i) SETERRQ(1,0,"IS is not sorted");
89 
90   ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr);
91   ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr);
92   ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr);
93   ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr);
94 
95   smap  = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap);
96   ssmap = smap;
97   lens  = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens);
98   PetscMemzero(smap,oldcols*sizeof(int));
99   for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1;
100   /* determine lens of each row */
101   for (i=0; i<nrows; i++) {
102     kstart  = ai[irow[i]];
103     kend    = kstart + a->ilen[irow[i]];
104     lens[i] = 0;
105       for ( k=kstart; k<kend; k++ ) {
106         if (ssmap[aj[k]]) {
107           lens[i]++;
108         }
109       }
110     }
111   /* Create and fill new matrix */
112   if (scall == MAT_REUSE_MATRIX) {
113     c = (Mat_SeqBAIJ *)((*B)->data);
114 
115     if (c->mbs!=nrows || c->nbs!=ncols || c->bs!=bs) SETERRQ(1,0,"");
116     if (PetscMemcmp(c->ilen,lens, c->mbs *sizeof(int))) {
117       SETERRQ(1,0,"Cannot reuse matrix. wrong no of nonzeros");
118     }
119     PetscMemzero(c->ilen,c->mbs*sizeof(int));
120     C = *B;
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   PetscFunctionReturn(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   PetscFunctionBegin;
164   ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr);
165   ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr);
166   ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr);
167   ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr);
168 
169   /* Verify if the indices corespond to each elementin a block
170    and form the IS with compressed IS */
171   vary = (int *) PetscMalloc(2*(a->mbs+1)*sizeof(int)); CHKPTRQ(vary);
172   iary = vary + a->mbs;
173   PetscMemzero(vary,(a->mbs)*sizeof(int));
174   for ( i=0; i<nrows; i++) vary[irow[i]/bs]++;
175   count = 0;
176   for (i=0; i<a->mbs; i++) {
177     if (vary[i]!=0 && vary[i]!=bs) SETERRA(1,0,"MatGetSubmatrices_SeqBAIJ:");
178     if (vary[i]==bs) iary[count++] = i;
179   }
180   ierr = ISCreateGeneral(PETSC_COMM_SELF, count, iary,&is1); CHKERRQ(ierr);
181 
182   PetscMemzero(vary,(a->mbs)*sizeof(int));
183   for ( i=0; i<ncols; i++) vary[icol[i]/bs]++;
184   count = 0;
185   for (i=0; i<a->mbs; i++) {
186     if (vary[i]!=0 && vary[i]!=bs) SETERRA(1,0,"MatGetSubmatrices_SeqBAIJ:");
187     if (vary[i]==bs) iary[count++] = i;
188   }
189   ierr = ISCreateGeneral(PETSC_COMM_SELF, count, iary,&is2); CHKERRQ(ierr);
190   ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr);
191   ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr);
192   PetscFree(vary);
193 
194   ierr = MatGetSubMatrix_SeqBAIJ_Private(A,is1,is2,scall,B); CHKERRQ(ierr);
195   ISDestroy(is1);
196   ISDestroy(is2);
197   PetscFunctionReturn(0);
198 }
199 
200 extern int MatGetSubMatrix_SeqBAIJ(Mat,IS,IS,MatGetSubMatrixCall,Mat*);
201 
202 #undef __FUNC__
203 #define __FUNC__ "MatGetSubMatrices_SeqBAIJ"
204 int MatGetSubMatrices_SeqBAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,
205                                     Mat **B)
206 {
207   int ierr,i;
208 
209   PetscFunctionBegin;
210   if (scall == MAT_INITIAL_MATRIX) {
211     *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B);
212   }
213 
214   for ( i=0; i<n; i++ ) {
215     ierr = MatGetSubMatrix_SeqBAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr);
216   }
217   PetscFunctionReturn(0);
218 }
219 
220 
221 
222 
223 
224 
225 
226