xref: /petsc/src/mat/impls/baij/seq/baijfact3.c (revision 098c0d8dc3b092b2dbccb3ec988c33e9401e6199) 
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] - bi[row] + 1; /* num of entries in the row with column index <= row */
174       nzbd    = bdiag[row] + 1; /* num of entries in the row with column index <= row */
175       ajtmp   = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */
176       ierr = PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);CHKERRQ(ierr);
177       nzi += nlnk;
178       row  = lnk[row];
179     }
180     bi[i+1] = bi[i] + nzi;
181     im[i]   = nzi;
182 
183     /* mark bdiag */
184     nzbd = 0;
185     nnz  = nzi;
186     k    = lnk[n];
187     while (nnz-- && k < i){
188       nzbd++;
189       k = lnk[k];
190     }
191     bdiag[i] = nzbd; /* note : bdaig[i] = nnzL as input for PetscFreeSpaceContiguous_newdatastruct() */
192 
193     /* if free space is not available, make more free space */
194     if (current_space->local_remaining<nzi) {
195       nnz = 2*(n - i)*nzi; /* estimated and max additional space needed */
196       ierr = PetscFreeSpaceGet(nnz,&current_space);CHKERRQ(ierr);
197       reallocs++;
198     }
199 
200     /* copy data into free space, then initialize lnk */
201     ierr = PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
202     bi_ptr[i] = current_space->array;
203     current_space->array           += nzi;
204     current_space->local_used      += nzi;
205     current_space->local_remaining -= nzi;
206   }
207 #if defined(PETSC_USE_INFO)
208   if (ai[n] != 0) {
209     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
210     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,f,af);CHKERRQ(ierr);
211     ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
212     ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G);\n",af);CHKERRQ(ierr);
213     ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr);
214   } else {
215     ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr);
216   }
217 #endif
218 
219   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
220   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
221 
222   /* destroy list of free space and other temporary array(s) */
223   ierr = PetscMalloc((bi[n]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
224   ierr = PetscFreeSpaceContiguous_newdatastruct(&free_space,bj,n,bi,bdiag);CHKERRQ(ierr);
225   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
226   ierr = PetscFree2(bi_ptr,im);CHKERRQ(ierr);
227 
228   /* put together the new matrix */
229   ierr = MatSeqBAIJSetPreallocation_SeqBAIJ(B,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr);
230   ierr = PetscLogObjectParent(B,isicol);CHKERRQ(ierr);
231   b    = (Mat_SeqBAIJ*)(B)->data;
232   b->free_a       = PETSC_TRUE;
233   b->free_ij      = PETSC_TRUE;
234   b->singlemalloc = PETSC_FALSE;
235   ierr          = PetscMalloc((bi[2*n+1])*sizeof(MatScalar)*bs2,&b->a);CHKERRQ(ierr);
236   b->j          = bj;
237   b->i          = bi;
238   b->diag       = bdiag;
239   b->ilen       = 0;
240   b->imax       = 0;
241   b->row        = isrow;
242   b->col        = iscol;
243   b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
244   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
245   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
246   b->icol       = isicol;
247   ierr = PetscMalloc((bs*n+bs)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
248   ierr = PetscLogObjectMemory(B,(bi[2*n+1])*(sizeof(PetscInt)+sizeof(PetscScalar)*bs2));CHKERRQ(ierr);
249 
250   b->maxnz = b->nz = bi[2*n+1] ;
251   (B)->factor                =  MAT_FACTOR_LU;
252   (B)->info.factor_mallocs   = reallocs;
253   (B)->info.fill_ratio_given = f;
254 
255   if (ai[n] != 0) {
256     (B)->info.fill_ratio_needed = ((PetscReal)bi[2*n+1])/((PetscReal)ai[n]);
257   } else {
258     (B)->info.fill_ratio_needed = 0.0;
259   }
260   (B)->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N_newdatastruct;
261 
262   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
263   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
264   both_identity = (PetscTruth) (row_identity && col_identity);
265   if(both_identity){
266      switch(bs){
267      case 2:
268         (B)->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering_newdatastruct;
269         break;
270      case 3:
271         (B)->ops->solve = MatSolve_SeqBAIJ_3_NaturalOrdering_newdatastruct;
272         break;
273      case 4:
274         (B)->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering_newdatastruct;
275         break;
276      case 5:
277         (B)->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering_newdatastruct;
278         break;
279      case 6:
280         (B)->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering_newdatastruct;
281         break;
282      case 7:
283         (B)->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering_newdatastruct;
284         break;
285      default:
286         (B)->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering_newdatastruct;
287       }
288    }
289    else {
290      switch(bs){
291      case 2:
292         (B)->ops->solve = MatSolve_SeqBAIJ_2_newdatastruct;
293         break;
294      case 3:
295         (B)->ops->solve = MatSolve_SeqBAIJ_3_newdatastruct;
296         break;
297      case 4:
298         (B)->ops->solve = MatSolve_SeqBAIJ_4_newdatastruct;
299         break;
300      case 5:
301         (B)->ops->solve = MatSolve_SeqBAIJ_5_newdatastruct;
302         break;
303      case 6:
304         (B)->ops->solve = MatSolve_SeqBAIJ_6_newdatastruct;
305         break;
306      case 7:
307         (B)->ops->solve = MatSolve_SeqBAIJ_7_newdatastruct;
308         break;
309      default:
310         (B)->ops->solve = MatSolve_SeqBAIJ_N_newdatastruct;
311       }
312    }
313 //  ierr = MatSeqBAIJSetNumericFactorization(B,both_identity);CHKERRQ(ierr);
314 
315   PetscFunctionReturn(0);
316  }
317 
318 #undef __FUNCT__
319 #define __FUNCT__ "MatLUFactorSymbolic_SeqBAIJ"
320 PetscErrorCode MatLUFactorSymbolic_SeqBAIJ(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
321 {
322   Mat_SeqBAIJ        *a = (Mat_SeqBAIJ*)A->data,*b;
323   PetscInt           n=a->mbs,bs = A->rmap->bs,bs2=a->bs2;
324   PetscTruth         row_identity,col_identity,both_identity;
325   IS                 isicol;
326   PetscErrorCode     ierr;
327   const PetscInt     *r,*ic;
328   PetscInt           i,*ai=a->i,*aj=a->j;
329   PetscInt           *bi,*bj,*ajtmp;
330   PetscInt           *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im;
331   PetscReal          f;
332   PetscInt           nlnk,*lnk,k,**bi_ptr;
333   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
334   PetscBT            lnkbt;
335   PetscTruth         newdatastruct=PETSC_FALSE;
336 
337   PetscFunctionBegin;
338   ierr = PetscOptionsGetTruth(PETSC_NULL,"-lu_new",&newdatastruct,PETSC_NULL);CHKERRQ(ierr);
339   if(newdatastruct){
340     ierr = MatLUFactorSymbolic_SeqBAIJ_newdatastruct(B,A,isrow,iscol,info);CHKERRQ(ierr);
341     PetscFunctionReturn(0);
342   }
343 
344   if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
345   ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr);
346   ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr);
347   ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr);
348 
349   /* get new row pointers */
350   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bi);CHKERRQ(ierr);
351   bi[0] = 0;
352 
353   /* bdiag is location of diagonal in factor */
354   ierr = PetscMalloc((n+1)*sizeof(PetscInt),&bdiag);CHKERRQ(ierr);
355   bdiag[0] = 0;
356 
357   /* linked list for storing column indices of the active row */
358   nlnk = n + 1;
359   ierr = PetscLLCreate(n,n,nlnk,lnk,lnkbt);CHKERRQ(ierr);
360 
361   ierr = PetscMalloc2(n+1,PetscInt**,&bi_ptr,n+1,PetscInt,&im);CHKERRQ(ierr);
362 
363   /* initial FreeSpace size is f*(ai[n]+1) */
364   f = info->fill;
365   ierr = PetscFreeSpaceGet((PetscInt)(f*(ai[n]+1)),&free_space);CHKERRQ(ierr);
366   current_space = free_space;
367 
368   for (i=0; i<n; i++) {
369     /* copy previous fill into linked list */
370     nzi = 0;
371     nnz = ai[r[i]+1] - ai[r[i]];
372     if (!nnz) SETERRQ2(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i);
373     ajtmp = aj + ai[r[i]];
374     ierr = PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);CHKERRQ(ierr);
375     nzi += nlnk;
376 
377     /* add pivot rows into linked list */
378     row = lnk[n];
379     while (row < i) {
380       nzbd    = bdiag[row] - bi[row] + 1; /* num of entries in the row with column index <= row */
381       ajtmp   = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */
382       ierr = PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);CHKERRQ(ierr);
383       nzi += nlnk;
384       row  = lnk[row];
385     }
386     bi[i+1] = bi[i] + nzi;
387     im[i]   = nzi;
388 
389     /* mark bdiag */
390     nzbd = 0;
391     nnz  = nzi;
392     k    = lnk[n];
393     while (nnz-- && k < i){
394       nzbd++;
395       k = lnk[k];
396     }
397     bdiag[i] = bi[i] + nzbd;
398 
399     /* if free space is not available, make more free space */
400     if (current_space->local_remaining<nzi) {
401       nnz = (n - i)*nzi; /* estimated and max additional space needed */
402       ierr = PetscFreeSpaceGet(nnz,&current_space);CHKERRQ(ierr);
403       reallocs++;
404     }
405 
406     /* copy data into free space, then initialize lnk */
407     ierr = PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
408     bi_ptr[i] = current_space->array;
409     current_space->array           += nzi;
410     current_space->local_used      += nzi;
411     current_space->local_remaining -= nzi;
412   }
413 #if defined(PETSC_USE_INFO)
414   if (ai[n] != 0) {
415     PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]);
416     ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %G needed %G\n",reallocs,f,af);CHKERRQ(ierr);
417     ierr = PetscInfo1(A,"Run with -pc_factor_fill %G or use \n",af);CHKERRQ(ierr);
418     ierr = PetscInfo1(A,"PCFactorSetFill(pc,%G);\n",af);CHKERRQ(ierr);
419     ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr);
420   } else {
421     ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr);
422   }
423 #endif
424 
425   ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr);
426   ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr);
427 
428   /* destroy list of free space and other temporary array(s) */
429   ierr = PetscMalloc((bi[n]+1)*sizeof(PetscInt),&bj);CHKERRQ(ierr);
430   ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr);
431   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
432   ierr = PetscFree2(bi_ptr,im);CHKERRQ(ierr);
433 
434   /* put together the new matrix */
435   ierr = MatSeqBAIJSetPreallocation_SeqBAIJ(B,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr);
436   ierr = PetscLogObjectParent(B,isicol);CHKERRQ(ierr);
437   b    = (Mat_SeqBAIJ*)(B)->data;
438   b->free_a       = PETSC_TRUE;
439   b->free_ij      = PETSC_TRUE;
440   b->singlemalloc = PETSC_FALSE;
441   ierr          = PetscMalloc((bi[n]+1)*sizeof(MatScalar)*bs2,&b->a);CHKERRQ(ierr);
442   b->j          = bj;
443   b->i          = bi;
444   b->diag       = bdiag;
445   b->ilen       = 0;
446   b->imax       = 0;
447   b->row        = isrow;
448   b->col        = iscol;
449   b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
450   ierr          = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr);
451   ierr          = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr);
452   b->icol       = isicol;
453   ierr = PetscMalloc((bs*n+bs)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
454   ierr = PetscLogObjectMemory(B,(bi[n]-n)*(sizeof(PetscInt)+sizeof(PetscScalar)*bs2));CHKERRQ(ierr);
455 
456   b->maxnz = b->nz = bi[n] ;
457   (B)->factor                =  MAT_FACTOR_LU;
458   (B)->info.factor_mallocs   = reallocs;
459   (B)->info.fill_ratio_given = f;
460 
461   if (ai[n] != 0) {
462     (B)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]);
463   } else {
464     (B)->info.fill_ratio_needed = 0.0;
465   }
466 
467   ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr);
468   ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr);
469   both_identity = (PetscTruth) (row_identity && col_identity);
470   ierr = MatSeqBAIJSetNumericFactorization(B,both_identity);CHKERRQ(ierr);
471   PetscFunctionReturn(0);
472  }
473 
474