xref: /petsc/src/mat/impls/baij/seq/baijfact3.c (revision b588c5a2f37da9e74d7d91cdaec80d3d0decb4d9) 
1 #define PETSCMAT_DLL
2 
3 /*
4     Factorization code for BAIJ format.
5 */
6 #include "../src/mat/impls/baij/seq/baij.h"
7 #include "../src/mat/blockinvert.h"
8 
9 #undef __FUNCT__
10 #define __FUNCT__ "MatSeqBAIJSetNumericFactorization"
11 /*
12    This is used to set the numeric factorization for both LU and ILU symbolic factorization
13 */
14 PetscErrorCode MatSeqBAIJSetNumericFactorization(Mat inA,PetscTruth natural)
15 {
16   PetscFunctionBegin;
17   if (natural) {
18     switch (inA->rmap->bs) {
19     case 1:
20       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1;
21       break;
22     case 2:
23       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering;
24       break;
25     case 3:
26       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering;
27       break;
28     case 4:
29 #if defined(PETSC_USE_MAT_SINGLE)
30       {
31         PetscTruth  sse_enabled_local;
32         PetscErrorCode ierr;
33         ierr = PetscSSEIsEnabled(inA->comm,&sse_enabled_local,PETSC_NULL);CHKERRQ(ierr);
34         if (sse_enabled_local) {
35 #  if defined(PETSC_HAVE_SSE)
36           int i,*AJ=a->j,nz=a->nz,n=a->mbs;
37           if (n==(unsigned short)n) {
38             unsigned short *aj=(unsigned short *)AJ;
39             for (i=0;i<nz;i++) {
40               aj[i] = (unsigned short)AJ[i];
41             }
42             inA->ops->setunfactored   = MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj;
43             inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering_SSE_usj;
44             ierr = PetscInfo(inA,"Using special SSE, in-place natural ordering, ushort j index factor BS=4\n");CHKERRQ(ierr);
45           } else {
46         /* Scale the column indices for easier indexing in MatSolve. */
47 /*            for (i=0;i<nz;i++) { */
48 /*              AJ[i] = AJ[i]*4; */
49 /*            } */
50             inA->ops->setunfactored   = MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE;
51             inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering_SSE;
52             ierr = PetscInfo(inA,"Using special SSE, in-place natural ordering, int j index factor BS=4\n");CHKERRQ(ierr);
53           }
54 #  else
55         /* This should never be reached.  If so, problem in PetscSSEIsEnabled. */
56           SETERRQ(PETSC_ERR_SUP,"SSE Hardware unavailable");
57 #  endif
58         } else {
59           inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering;
60         }
61       }
62 #else
63       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering;
64 #endif
65       break;
66     case 5:
67       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering;
68       break;
69     case 6:
70       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering;
71       break;
72     case 7:
73       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering;
74       break;
75     default:
76       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N;
77       break;
78     }
79   } else {
80     switch (inA->rmap->bs) {
81     case 1:
82       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1;
83       break;
84     case 2:
85       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2;
86       break;
87     case 3:
88       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3;
89       break;
90     case 4:
91       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4;
92       break;
93     case 5:
94       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5;
95       break;
96     case 6:
97       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6;
98       break;
99     case 7:
100       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7;
101       break;
102     default:
103       inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N;
104       break;
105     }
106   }
107   PetscFunctionReturn(0);
108 }
109 
110 /*
111     The symbolic factorization code is identical to that for AIJ format,
112   except for very small changes since this is now a SeqBAIJ datastructure.
113   NOT good code reuse.
114 */
115 #include "petscbt.h"
116 #include "../src/mat/utils/freespace.h"
117 
118 #undef __FUNCT__
119 #define __FUNCT__ "MatLUFactorSymbolic_SeqBAIJ_newdatastruct"
120 PetscErrorCode MatLUFactorSymbolic_SeqBAIJ_newdatastruct(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
121 {
122   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data,*b;
123   PetscInt           n=a->mbs,bs = A->rmap->bs,bs2=a->bs2;
124   PetscTruth         row_identity,col_identity,both_identity;
125   IS                 isicol;
126   PetscErrorCode     ierr;
127   const PetscInt     *r,*ic;
128   PetscInt           i,*ai=a->i,*aj=a->j;
129   PetscInt           *bi,*bj,*ajtmp;
130   PetscInt           *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im;
131   PetscReal          f;
132   PetscInt           nlnk,*lnk,k,**bi_ptr;
133   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
134   PetscBT            lnkbt;
135 
136   PetscFunctionBegin;
137   if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
138   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
139   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
140   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
141 
142   /* get new row pointers */
143   ierr = PetscMalloc((2*n+2)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
144   bi[0] = 0;
145 
146   /* bdiag is location of diagonal in factor */
147   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bdiag);CHKERRQ(ierr);
148   bdiag[0] = 0;
149 
150   /* linked list for storing column indices of the active row */
151   nlnk = n + 1;
152   ierr = PetscLLCreate(n,n,nlnk,lnk,lnkbt);CHKERRQ(ierr);
153 
154   ierr = PetscMalloc2(n+1,PetscInt**,&bi_ptr,n+1,PetscInt,&im);CHKERRQ(ierr);
155 
156   /* initial FreeSpace size is f*(ai[n]+1) */
157   f = info->fill;
158   ierr = PetscFreeSpaceGet((PetscInt)(f*(ai[n]+1)),&free_space);CHKERRQ(ierr);
159   current_space = free_space;
160 
161   for (i=0; i<n; i++) {
162     /* copy previous fill into linked list */
163     nzi = 0;
164     nnz = ai[r[i]+1] - ai[r[i]];
165     if (!nnz) SETERRQ2(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i);
166     ajtmp = aj + ai[r[i]];
167     ierr = PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);CHKERRQ(ierr);
168     nzi += nlnk;
169 
170     /* add pivot rows into linked list */
171     row = lnk[n];
172     while (row < i) {
173       nzbd    = bdiag[row] + 1; /* num of entries in the row with column index <= row */
174       ajtmp   = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */
175       ierr = PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);CHKERRQ(ierr);
176       nzi += nlnk;
177       row  = lnk[row];
178     }
179     bi[i+1] = bi[i] + nzi;
180     im[i]   = nzi;
181 
182     /* mark bdiag */
183     nzbd = 0;
184     nnz  = nzi;
185     k    = lnk[n];
186     while (nnz-- && k < i){
187       nzbd++;
188       k = lnk[k];
189     }
190     bdiag[i] = nzbd; /* note : bdaig[i] = nnzL as input for PetscFreeSpaceContiguous_LU() */
191 
192     /* if free space is not available, make more free space */
193     if (current_space->local_remaining<nzi) {
194       nnz = 2*(n - i)*nzi; /* estimated and max additional space needed */
195       ierr = PetscFreeSpaceGet(nnz,&current_space);CHKERRQ(ierr);
196       reallocs++;
197     }
198 
199     /* copy data into free space, then initialize lnk */
200     ierr = PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
201     bi_ptr[i] = current_space->array;
202     current_space->array           += nzi;
203     current_space->local_used      += nzi;
204     current_space->local_remaining -= nzi;
205   }
206 #if defined(PETSC_USE_INFO)
207   if (ai[n] != 0) {
208     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
209     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,f,af);CHKERRQ(ierr);
210     ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
211     ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G);\n",af);CHKERRQ(ierr);
212     ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr);
213   } else {
214     ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr);
215   }
216 #endif
217 
218   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
219   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
220 
221   /* destroy list of free space and other temporary array(s) */
222   ierr = PetscMalloc((bi[n]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
223   ierr = PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);CHKERRQ(ierr);
224   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
225   ierr = PetscFree2(bi_ptr,im);CHKERRQ(ierr);
226 
227   /* put together the new matrix */
228   ierr = MatSeqBAIJSetPreallocation_SeqBAIJ(B,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr);
229   ierr = PetscLogObjectParent(B,isicol);CHKERRQ(ierr);
230   b    = (Mat_SeqBAIJ*)(B)->data;
231   b->free_a       = PETSC_TRUE;
232   b->free_ij      = PETSC_TRUE;
233   b->singlemalloc = PETSC_FALSE;
234   ierr          = PetscMalloc((bi[2*n+1])*sizeof(MatScalar)*bs2,&b->a);CHKERRQ(ierr);
235   b->j          = bj;
236   b->i          = bi;
237   b->diag       = bdiag;
238   b->ilen       = 0;
239   b->imax       = 0;
240   b->row        = isrow;
241   b->col        = iscol;
242   b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
243   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
244   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
245   b->icol       = isicol;
246   ierr = PetscMalloc((bs*n+bs)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
247   ierr = PetscLogObjectMemory(B,(bi[2*n+1])*(sizeof(PetscInt)+sizeof(PetscScalar)*bs2));CHKERRQ(ierr);
248 
249   b->maxnz = b->nz = bi[2*n+1];
250   (B)->factor                =  MAT_FACTOR_LU;
251   (B)->info.factor_mallocs   = reallocs;
252   (B)->info.fill_ratio_given = f;
253 
254   if (ai[n] != 0) {
255     (B)->info.fill_ratio_needed = ((PetscReal)bi[2*n+1])/((PetscReal)ai[n]);
256   } else {
257     (B)->info.fill_ratio_needed = 0.0;
258   }
259 
260   switch(bs){
261     case 2:
262       (B)->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_newdatastruct;
263       break;
264     default:
265       (B)->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N_newdatastruct;
266   }
267 
268   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
269   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
270   both_identity = (PetscTruth) (row_identity && col_identity);
271   if(both_identity){
272      switch(bs){
273      case 2:
274         (B)->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering_newdatastruct;
275         break;
276      case 3:
277         (B)->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering_newdatastruct;
278         break;
279      case 4:
280         (B)->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering_newdatastruct;
281         break;
282      case 5:
283         (B)->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering_newdatastruct;
284         break;
285      case 6:
286         (B)->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering_newdatastruct;
287         break;
288      case 7:
289         (B)->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering_newdatastruct;
290         break;
291      default:
292         (B)->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering_newdatastruct;
293       }
294    }
295    else {
296      switch(bs){
297      case 2:
298         (B)->ops->solve = MatSolve_SeqBAIJ_2_newdatastruct;
299         break;
300      case 3:
301         (B)->ops->solve = MatSolve_SeqBAIJ_3_newdatastruct;
302         break;
303      case 4:
304         (B)->ops->solve = MatSolve_SeqBAIJ_4_newdatastruct;
305         break;
306      case 5:
307         (B)->ops->solve = MatSolve_SeqBAIJ_5_newdatastruct;
308         break;
309      case 6:
310         (B)->ops->solve = MatSolve_SeqBAIJ_6_newdatastruct;
311         break;
312      case 7:
313         (B)->ops->solve = MatSolve_SeqBAIJ_7_newdatastruct;
314         break;
315      default:
316         (B)->ops->solve = MatSolve_SeqBAIJ_N_newdatastruct;
317       }
318    }
319   /*  ierr = MatSeqBAIJSetNumericFactorization(B,both_identity);CHKERRQ(ierr); */
320   PetscFunctionReturn(0);
321  }
322 
323 #undef __FUNCT__
324 #define __FUNCT__ "MatLUFactorSymbolic_SeqBAIJ"
325 PetscErrorCode MatLUFactorSymbolic_SeqBAIJ(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
326 {
327   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data,*b;
328   PetscInt           n=a->mbs,bs = A->rmap->bs,bs2=a->bs2;
329   PetscTruth         row_identity,col_identity,both_identity;
330   IS                 isicol;
331   PetscErrorCode     ierr;
332   const PetscInt     *r,*ic;
333   PetscInt           i,*ai=a->i,*aj=a->j;
334   PetscInt           *bi,*bj,*ajtmp;
335   PetscInt           *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im;
336   PetscReal          f;
337   PetscInt           nlnk,*lnk,k,**bi_ptr;
338   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
339   PetscBT            lnkbt;
340   PetscTruth         newdatastruct=PETSC_FALSE;
341 
342   PetscFunctionBegin;
343   ierr = PetscOptionsGetTruth(PETSC_NULL,"-lu_new",&newdatastruct,PETSC_NULL);CHKERRQ(ierr);
344   if(newdatastruct){
345     ierr = MatLUFactorSymbolic_SeqBAIJ_newdatastruct(B,A,isrow,iscol,info);CHKERRQ(ierr);
346     PetscFunctionReturn(0);
347   }
348 
349   if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
350   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
351   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
352   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
353 
354   /* get new row pointers */
355   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
356   bi[0] = 0;
357 
358   /* bdiag is location of diagonal in factor */
359   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bdiag);CHKERRQ(ierr);
360   bdiag[0] = 0;
361 
362   /* linked list for storing column indices of the active row */
363   nlnk = n + 1;
364   ierr = PetscLLCreate(n,n,nlnk,lnk,lnkbt);CHKERRQ(ierr);
365 
366   ierr = PetscMalloc2(n+1,PetscInt**,&bi_ptr,n+1,PetscInt,&im);CHKERRQ(ierr);
367 
368   /* initial FreeSpace size is f*(ai[n]+1) */
369   f = info->fill;
370   ierr = PetscFreeSpaceGet((PetscInt)(f*(ai[n]+1)),&free_space);CHKERRQ(ierr);
371   current_space = free_space;
372 
373   for (i=0; i<n; i++) {
374     /* copy previous fill into linked list */
375     nzi = 0;
376     nnz = ai[r[i]+1] - ai[r[i]];
377     if (!nnz) SETERRQ2(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i);
378     ajtmp = aj + ai[r[i]];
379     ierr = PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);CHKERRQ(ierr);
380     nzi += nlnk;
381 
382     /* add pivot rows into linked list */
383     row = lnk[n];
384     while (row < i) {
385       nzbd    = bdiag[row] - bi[row] + 1; /* num of entries in the row with column index <= row */
386       ajtmp   = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */
387       ierr = PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);CHKERRQ(ierr);
388       nzi += nlnk;
389       row  = lnk[row];
390     }
391     bi[i+1] = bi[i] + nzi;
392     im[i]   = nzi;
393 
394     /* mark bdiag */
395     nzbd = 0;
396     nnz  = nzi;
397     k    = lnk[n];
398     while (nnz-- && k < i){
399       nzbd++;
400       k = lnk[k];
401     }
402     bdiag[i] = bi[i] + nzbd;
403 
404     /* if free space is not available, make more free space */
405     if (current_space->local_remaining<nzi) {
406       nnz = (n - i)*nzi; /* estimated and max additional space needed */
407       ierr = PetscFreeSpaceGet(nnz,&current_space);CHKERRQ(ierr);
408       reallocs++;
409     }
410 
411     /* copy data into free space, then initialize lnk */
412     ierr = PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
413     bi_ptr[i] = current_space->array;
414     current_space->array           += nzi;
415     current_space->local_used      += nzi;
416     current_space->local_remaining -= nzi;
417   }
418 #if defined(PETSC_USE_INFO)
419   if (ai[n] != 0) {
420     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
421     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,f,af);CHKERRQ(ierr);
422     ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
423     ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G);\n",af);CHKERRQ(ierr);
424     ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr);
425   } else {
426     ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr);
427   }
428 #endif
429 
430   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
431   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
432 
433   /* destroy list of free space and other temporary array(s) */
434   ierr = PetscMalloc((bi[n]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
435   ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr);
436   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
437   ierr = PetscFree2(bi_ptr,im);CHKERRQ(ierr);
438 
439   /* put together the new matrix */
440   ierr = MatSeqBAIJSetPreallocation_SeqBAIJ(B,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr);
441   ierr = PetscLogObjectParent(B,isicol);CHKERRQ(ierr);
442   b    = (Mat_SeqBAIJ*)(B)->data;
443   b->free_a       = PETSC_TRUE;
444   b->free_ij      = PETSC_TRUE;
445   b->singlemalloc = PETSC_FALSE;
446   ierr          = PetscMalloc((bi[n]+1)*sizeof(MatScalar)*bs2,&b->a);CHKERRQ(ierr);
447   b->j          = bj;
448   b->i          = bi;
449   b->diag       = bdiag;
450   b->ilen       = 0;
451   b->imax       = 0;
452   b->row        = isrow;
453   b->col        = iscol;
454   b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
455   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
456   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
457   b->icol       = isicol;
458   ierr = PetscMalloc((bs*n+bs)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
459   ierr = PetscLogObjectMemory(B,(bi[n]-n)*(sizeof(PetscInt)+sizeof(PetscScalar)*bs2));CHKERRQ(ierr);
460 
461   b->maxnz = b->nz = bi[n] ;
462   (B)->factor                =  MAT_FACTOR_LU;
463   (B)->info.factor_mallocs   = reallocs;
464   (B)->info.fill_ratio_given = f;
465 
466   if (ai[n] != 0) {
467     (B)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
468   } else {
469     (B)->info.fill_ratio_needed = 0.0;
470   }
471 
472   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
473   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
474   both_identity = (PetscTruth) (row_identity && col_identity);
475   ierr = MatSeqBAIJSetNumericFactorization(B,both_identity);CHKERRQ(ierr);
476   PetscFunctionReturn(0);
477  }
478 
479