Actual source code: baij.c

petsc-3.13.6 2020-09-29
Report Typos and Errors

  2: /*
  3:     Defines the basic matrix operations for the BAIJ (compressed row)
  4:   matrix storage format.
  5: */
  6:  #include <../src/mat/impls/baij/seq/baij.h>
  7:  #include <petscblaslapack.h>
  8:  #include <petsc/private/kernels/blockinvert.h>
  9:  #include <petsc/private/kernels/blockmatmult.h>

 11: #if defined(PETSC_HAVE_HYPRE)
 12: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat,MatType,MatReuse,Mat*);
 13: #endif

 15: #if defined(PETSC_HAVE_MKL_SPARSE_OPTIMIZE)
 16: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqBAIJMKL(Mat,MatType,MatReuse,Mat*);
 17: #endif
 18: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*);

 20: PetscErrorCode MatInvertBlockDiagonal_SeqBAIJ(Mat A,const PetscScalar **values)
 21: {
 22:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*) A->data;
 24:   PetscInt       *diag_offset,i,bs = A->rmap->bs,mbs = a->mbs,ipvt[5],bs2 = bs*bs,*v_pivots;
 25:   MatScalar      *v    = a->a,*odiag,*diag,work[25],*v_work;
 26:   PetscReal      shift = 0.0;
 27:   PetscBool      allowzeropivot,zeropivotdetected=PETSC_FALSE;

 30:   allowzeropivot = PetscNot(A->erroriffailure);

 32:   if (a->idiagvalid) {
 33:     if (values) *values = a->idiag;
 34:     return(0);
 35:   }
 36:   MatMarkDiagonal_SeqBAIJ(A);
 37:   diag_offset = a->diag;
 38:   if (!a->idiag) {
 39:     PetscMalloc1(bs2*mbs,&a->idiag);
 40:     PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));
 41:   }
 42:   diag  = a->idiag;
 43:   if (values) *values = a->idiag;
 44:   /* factor and invert each block */
 45:   switch (bs) {
 46:   case 1:
 47:     for (i=0; i<mbs; i++) {
 48:       odiag    = v + 1*diag_offset[i];
 49:       diag[0]  = odiag[0];

 51:       if (PetscAbsScalar(diag[0] + shift) < PETSC_MACHINE_EPSILON) {
 52:         if (allowzeropivot) {
 53:           A->factorerrortype             = MAT_FACTOR_NUMERIC_ZEROPIVOT;
 54:           A->factorerror_zeropivot_value = PetscAbsScalar(diag[0]);
 55:           A->factorerror_zeropivot_row   = i;
 56:           PetscInfo1(A,"Zero pivot, row %D\n",i);
 57:         } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %D pivot value %g tolerance %g",i,(double)PetscAbsScalar(diag[0]),(double)PETSC_MACHINE_EPSILON);
 58:       }

 60:       diag[0]  = (PetscScalar)1.0 / (diag[0] + shift);
 61:       diag    += 1;
 62:     }
 63:     break;
 64:   case 2:
 65:     for (i=0; i<mbs; i++) {
 66:       odiag    = v + 4*diag_offset[i];
 67:       diag[0]  = odiag[0]; diag[1] = odiag[1]; diag[2] = odiag[2]; diag[3] = odiag[3];
 68:       PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);
 69:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
 70:       diag    += 4;
 71:     }
 72:     break;
 73:   case 3:
 74:     for (i=0; i<mbs; i++) {
 75:       odiag    = v + 9*diag_offset[i];
 76:       diag[0]  = odiag[0]; diag[1] = odiag[1]; diag[2] = odiag[2]; diag[3] = odiag[3];
 77:       diag[4]  = odiag[4]; diag[5] = odiag[5]; diag[6] = odiag[6]; diag[7] = odiag[7];
 78:       diag[8]  = odiag[8];
 79:       PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);
 80:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
 81:       diag    += 9;
 82:     }
 83:     break;
 84:   case 4:
 85:     for (i=0; i<mbs; i++) {
 86:       odiag  = v + 16*diag_offset[i];
 87:       PetscArraycpy(diag,odiag,16);
 88:       PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);
 89:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
 90:       diag  += 16;
 91:     }
 92:     break;
 93:   case 5:
 94:     for (i=0; i<mbs; i++) {
 95:       odiag  = v + 25*diag_offset[i];
 96:       PetscArraycpy(diag,odiag,25);
 97:       PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);
 98:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
 99:       diag  += 25;
100:     }
101:     break;
102:   case 6:
103:     for (i=0; i<mbs; i++) {
104:       odiag  = v + 36*diag_offset[i];
105:       PetscArraycpy(diag,odiag,36);
106:       PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);
107:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
108:       diag  += 36;
109:     }
110:     break;
111:   case 7:
112:     for (i=0; i<mbs; i++) {
113:       odiag  = v + 49*diag_offset[i];
114:       PetscArraycpy(diag,odiag,49);
115:       PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);
116:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
117:       diag  += 49;
118:     }
119:     break;
120:   default:
121:     PetscMalloc2(bs,&v_work,bs,&v_pivots);
122:     for (i=0; i<mbs; i++) {
123:       odiag  = v + bs2*diag_offset[i];
124:       PetscArraycpy(diag,odiag,bs2);
125:       PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);
126:       if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
127:       diag  += bs2;
128:     }
129:     PetscFree2(v_work,v_pivots);
130:   }
131:   a->idiagvalid = PETSC_TRUE;
132:   return(0);
133: }

135: PetscErrorCode MatSOR_SeqBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
136: {
137:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
138:   PetscScalar       *x,*work,*w,*workt,*t;
139:   const MatScalar   *v,*aa = a->a, *idiag;
140:   const PetscScalar *b,*xb;
141:   PetscScalar       s[7], xw[7]={0}; /* avoid some compilers thinking xw is uninitialized */
142:   PetscErrorCode    ierr;
143:   PetscInt          m = a->mbs,i,i2,nz,bs = A->rmap->bs,bs2 = bs*bs,k,j,idx,it;
144:   const PetscInt    *diag,*ai = a->i,*aj = a->j,*vi;

147:   its = its*lits;
148:   if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat");
149:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
150:   if (fshift) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Sorry, no support for diagonal shift");
151:   if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Sorry, no support for non-trivial relaxation factor");
152:   if ((flag & SOR_APPLY_UPPER) || (flag & SOR_APPLY_LOWER)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Sorry, no support for applying upper or lower triangular parts");

154:   if (!a->idiagvalid) {MatInvertBlockDiagonal(A,NULL);}

156:   if (!m) return(0);
157:   diag  = a->diag;
158:   idiag = a->idiag;
159:   k    = PetscMax(A->rmap->n,A->cmap->n);
160:   if (!a->mult_work) {
161:     PetscMalloc1(k+1,&a->mult_work);
162:   }
163:   if (!a->sor_workt) {
164:     PetscMalloc1(k,&a->sor_workt);
165:   }
166:   if (!a->sor_work) {
167:     PetscMalloc1(bs,&a->sor_work);
168:   }
169:   work = a->mult_work;
170:   t    = a->sor_workt;
171:   w    = a->sor_work;

173:   VecGetArray(xx,&x);
174:   VecGetArrayRead(bb,&b);

176:   if (flag & SOR_ZERO_INITIAL_GUESS) {
177:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
178:       switch (bs) {
179:       case 1:
180:         PetscKernel_v_gets_A_times_w_1(x,idiag,b);
181:         t[0] = b[0];
182:         i2     = 1;
183:         idiag += 1;
184:         for (i=1; i<m; i++) {
185:           v  = aa + ai[i];
186:           vi = aj + ai[i];
187:           nz = diag[i] - ai[i];
188:           s[0] = b[i2];
189:           for (j=0; j<nz; j++) {
190:             xw[0] = x[vi[j]];
191:             PetscKernel_v_gets_v_minus_A_times_w_1(s,(v+j),xw);
192:           }
193:           t[i2] = s[0];
194:           PetscKernel_v_gets_A_times_w_1(xw,idiag,s);
195:           x[i2]  = xw[0];
196:           idiag += 1;
197:           i2    += 1;
198:         }
199:         break;
200:       case 2:
201:         PetscKernel_v_gets_A_times_w_2(x,idiag,b);
202:         t[0] = b[0]; t[1] = b[1];
203:         i2     = 2;
204:         idiag += 4;
205:         for (i=1; i<m; i++) {
206:           v  = aa + 4*ai[i];
207:           vi = aj + ai[i];
208:           nz = diag[i] - ai[i];
209:           s[0] = b[i2]; s[1] = b[i2+1];
210:           for (j=0; j<nz; j++) {
211:             idx = 2*vi[j];
212:             it  = 4*j;
213:             xw[0] = x[idx]; xw[1] = x[1+idx];
214:             PetscKernel_v_gets_v_minus_A_times_w_2(s,(v+it),xw);
215:           }
216:           t[i2] = s[0]; t[i2+1] = s[1];
217:           PetscKernel_v_gets_A_times_w_2(xw,idiag,s);
218:           x[i2]   = xw[0]; x[i2+1] = xw[1];
219:           idiag  += 4;
220:           i2     += 2;
221:         }
222:         break;
223:       case 3:
224:         PetscKernel_v_gets_A_times_w_3(x,idiag,b);
225:         t[0] = b[0]; t[1] = b[1]; t[2] = b[2];
226:         i2     = 3;
227:         idiag += 9;
228:         for (i=1; i<m; i++) {
229:           v  = aa + 9*ai[i];
230:           vi = aj + ai[i];
231:           nz = diag[i] - ai[i];
232:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2];
233:           while (nz--) {
234:             idx = 3*(*vi++);
235:             xw[0] = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx];
236:             PetscKernel_v_gets_v_minus_A_times_w_3(s,v,xw);
237:             v  += 9;
238:           }
239:           t[i2] = s[0]; t[i2+1] = s[1]; t[i2+2] = s[2];
240:           PetscKernel_v_gets_A_times_w_3(xw,idiag,s);
241:           x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2];
242:           idiag  += 9;
243:           i2     += 3;
244:         }
245:         break;
246:       case 4:
247:         PetscKernel_v_gets_A_times_w_4(x,idiag,b);
248:         t[0] = b[0]; t[1] = b[1]; t[2] = b[2]; t[3] = b[3];
249:         i2     = 4;
250:         idiag += 16;
251:         for (i=1; i<m; i++) {
252:           v  = aa + 16*ai[i];
253:           vi = aj + ai[i];
254:           nz = diag[i] - ai[i];
255:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3];
256:           while (nz--) {
257:             idx = 4*(*vi++);
258:             xw[0]  = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx]; xw[3] = x[3+idx];
259:             PetscKernel_v_gets_v_minus_A_times_w_4(s,v,xw);
260:             v  += 16;
261:           }
262:           t[i2] = s[0]; t[i2+1] = s[1]; t[i2+2] = s[2]; t[i2 + 3] = s[3];
263:           PetscKernel_v_gets_A_times_w_4(xw,idiag,s);
264:           x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3];
265:           idiag  += 16;
266:           i2     += 4;
267:         }
268:         break;
269:       case 5:
270:         PetscKernel_v_gets_A_times_w_5(x,idiag,b);
271:         t[0] = b[0]; t[1] = b[1]; t[2] = b[2]; t[3] = b[3]; t[4] = b[4];
272:         i2     = 5;
273:         idiag += 25;
274:         for (i=1; i<m; i++) {
275:           v  = aa + 25*ai[i];
276:           vi = aj + ai[i];
277:           nz = diag[i] - ai[i];
278:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3]; s[4] = b[i2+4];
279:           while (nz--) {
280:             idx = 5*(*vi++);
281:             xw[0]  = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx]; xw[3] = x[3+idx]; xw[4] = x[4+idx];
282:             PetscKernel_v_gets_v_minus_A_times_w_5(s,v,xw);
283:             v  += 25;
284:           }
285:           t[i2] = s[0]; t[i2+1] = s[1]; t[i2+2] = s[2]; t[i2+3] = s[3]; t[i2+4] = s[4];
286:           PetscKernel_v_gets_A_times_w_5(xw,idiag,s);
287:           x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3]; x[i2+4] = xw[4];
288:           idiag  += 25;
289:           i2     += 5;
290:         }
291:         break;
292:       case 6:
293:         PetscKernel_v_gets_A_times_w_6(x,idiag,b);
294:         t[0] = b[0]; t[1] = b[1]; t[2] = b[2]; t[3] = b[3]; t[4] = b[4]; t[5] = b[5];
295:         i2     = 6;
296:         idiag += 36;
297:         for (i=1; i<m; i++) {
298:           v  = aa + 36*ai[i];
299:           vi = aj + ai[i];
300:           nz = diag[i] - ai[i];
301:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3]; s[4] = b[i2+4]; s[5] = b[i2+5];
302:           while (nz--) {
303:             idx = 6*(*vi++);
304:             xw[0] = x[idx];   xw[1] = x[1+idx]; xw[2] = x[2+idx];
305:             xw[3] = x[3+idx]; xw[4] = x[4+idx]; xw[5] = x[5+idx];
306:             PetscKernel_v_gets_v_minus_A_times_w_6(s,v,xw);
307:             v  += 36;
308:           }
309:           t[i2]   = s[0]; t[i2+1] = s[1]; t[i2+2] = s[2];
310:           t[i2+3] = s[3]; t[i2+4] = s[4]; t[i2+5] = s[5];
311:           PetscKernel_v_gets_A_times_w_6(xw,idiag,s);
312:           x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3]; x[i2+4] = xw[4]; x[i2+5] = xw[5];
313:           idiag  += 36;
314:           i2     += 6;
315:         }
316:         break;
317:       case 7:
318:         PetscKernel_v_gets_A_times_w_7(x,idiag,b);
319:         t[0] = b[0]; t[1] = b[1]; t[2] = b[2];
320:         t[3] = b[3]; t[4] = b[4]; t[5] = b[5]; t[6] = b[6];
321:         i2     = 7;
322:         idiag += 49;
323:         for (i=1; i<m; i++) {
324:           v  = aa + 49*ai[i];
325:           vi = aj + ai[i];
326:           nz = diag[i] - ai[i];
327:           s[0] = b[i2];   s[1] = b[i2+1]; s[2] = b[i2+2];
328:           s[3] = b[i2+3]; s[4] = b[i2+4]; s[5] = b[i2+5]; s[6] = b[i2+6];
329:           while (nz--) {
330:             idx = 7*(*vi++);
331:             xw[0] = x[idx];   xw[1] = x[1+idx]; xw[2] = x[2+idx];
332:             xw[3] = x[3+idx]; xw[4] = x[4+idx]; xw[5] = x[5+idx]; xw[6] = x[6+idx];
333:             PetscKernel_v_gets_v_minus_A_times_w_7(s,v,xw);
334:             v  += 49;
335:           }
336:           t[i2]   = s[0]; t[i2+1] = s[1]; t[i2+2] = s[2];
337:           t[i2+3] = s[3]; t[i2+4] = s[4]; t[i2+5] = s[5]; t[i2+6] = s[6];
338:           PetscKernel_v_gets_A_times_w_7(xw,idiag,s);
339:           x[i2] =   xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2];
340:           x[i2+3] = xw[3]; x[i2+4] = xw[4]; x[i2+5] = xw[5]; x[i2+6] = xw[6];
341:           idiag  += 49;
342:           i2     += 7;
343:         }
344:         break;
345:       default:
346:         PetscKernel_w_gets_Ar_times_v(bs,bs,b,idiag,x);
347:         PetscArraycpy(t,b,bs);
348:         i2     = bs;
349:         idiag += bs2;
350:         for (i=1; i<m; i++) {
351:           v  = aa + bs2*ai[i];
352:           vi = aj + ai[i];
353:           nz = diag[i] - ai[i];

355:           PetscArraycpy(w,b+i2,bs);
356:           /* copy all rows of x that are needed into contiguous space */
357:           workt = work;
358:           for (j=0; j<nz; j++) {
359:             PetscArraycpy(workt,x + bs*(*vi++),bs);
360:             workt += bs;
361:           }
362:           PetscKernel_w_gets_w_minus_Ar_times_v(bs,bs*nz,w,v,work);
363:           PetscArraycpy(t+i2,w,bs);
364:           PetscKernel_w_gets_Ar_times_v(bs,bs,w,idiag,x+i2);

366:           idiag += bs2;
367:           i2    += bs;
368:         }
369:         break;
370:       }
371:       /* for logging purposes assume number of nonzero in lower half is 1/2 of total */
372:       PetscLogFlops(1.0*bs2*a->nz);
373:       xb = t;
374:     }
375:     else xb = b;
376:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
377:       idiag = a->idiag+bs2*(a->mbs-1);
378:       i2 = bs * (m-1);
379:       switch (bs) {
380:       case 1:
381:         s[0]  = xb[i2];
382:         PetscKernel_v_gets_A_times_w_1(xw,idiag,s);
383:         x[i2] = xw[0];
384:         i2   -= 1;
385:         for (i=m-2; i>=0; i--) {
386:           v  = aa + (diag[i]+1);
387:           vi = aj + diag[i] + 1;
388:           nz = ai[i+1] - diag[i] - 1;
389:           s[0] = xb[i2];
390:           for (j=0; j<nz; j++) {
391:             xw[0] = x[vi[j]];
392:             PetscKernel_v_gets_v_minus_A_times_w_1(s,(v+j),xw);
393:           }
394:           PetscKernel_v_gets_A_times_w_1(xw,idiag,s);
395:           x[i2]  = xw[0];
396:           idiag -= 1;
397:           i2    -= 1;
398:         }
399:         break;
400:       case 2:
401:         s[0]  = xb[i2]; s[1] = xb[i2+1];
402:         PetscKernel_v_gets_A_times_w_2(xw,idiag,s);
403:         x[i2] = xw[0]; x[i2+1] = xw[1];
404:         i2    -= 2;
405:         idiag -= 4;
406:         for (i=m-2; i>=0; i--) {
407:           v  = aa + 4*(diag[i] + 1);
408:           vi = aj + diag[i] + 1;
409:           nz = ai[i+1] - diag[i] - 1;
410:           s[0] = xb[i2]; s[1] = xb[i2+1];
411:           for (j=0; j<nz; j++) {
412:             idx = 2*vi[j];
413:             it  = 4*j;
414:             xw[0] = x[idx]; xw[1] = x[1+idx];
415:             PetscKernel_v_gets_v_minus_A_times_w_2(s,(v+it),xw);
416:           }
417:           PetscKernel_v_gets_A_times_w_2(xw,idiag,s);
418:           x[i2]   = xw[0]; x[i2+1] = xw[1];
419:           idiag  -= 4;
420:           i2     -= 2;
421:         }
422:         break;
423:       case 3:
424:         s[0]  = xb[i2]; s[1] = xb[i2+1]; s[2] = xb[i2+2];
425:         PetscKernel_v_gets_A_times_w_3(xw,idiag,s);
426:         x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2];
427:         i2    -= 3;
428:         idiag -= 9;
429:         for (i=m-2; i>=0; i--) {
430:           v  = aa + 9*(diag[i]+1);
431:           vi = aj + diag[i] + 1;
432:           nz = ai[i+1] - diag[i] - 1;
433:           s[0] = xb[i2]; s[1] = xb[i2+1]; s[2] = xb[i2+2];
434:           while (nz--) {
435:             idx = 3*(*vi++);
436:             xw[0] = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx];
437:             PetscKernel_v_gets_v_minus_A_times_w_3(s,v,xw);
438:             v  += 9;
439:           }
440:           PetscKernel_v_gets_A_times_w_3(xw,idiag,s);
441:           x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2];
442:           idiag  -= 9;
443:           i2     -= 3;
444:         }
445:         break;
446:       case 4:
447:         s[0]  = xb[i2]; s[1] = xb[i2+1]; s[2] = xb[i2+2]; s[3] = xb[i2+3];
448:         PetscKernel_v_gets_A_times_w_4(xw,idiag,s);
449:         x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3];
450:         i2    -= 4;
451:         idiag -= 16;
452:         for (i=m-2; i>=0; i--) {
453:           v  = aa + 16*(diag[i]+1);
454:           vi = aj + diag[i] + 1;
455:           nz = ai[i+1] - diag[i] - 1;
456:           s[0] = xb[i2]; s[1] = xb[i2+1]; s[2] = xb[i2+2]; s[3] = xb[i2+3];
457:           while (nz--) {
458:             idx = 4*(*vi++);
459:             xw[0]  = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx]; xw[3] = x[3+idx];
460:             PetscKernel_v_gets_v_minus_A_times_w_4(s,v,xw);
461:             v  += 16;
462:           }
463:           PetscKernel_v_gets_A_times_w_4(xw,idiag,s);
464:           x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3];
465:           idiag  -= 16;
466:           i2     -= 4;
467:         }
468:         break;
469:       case 5:
470:         s[0]  = xb[i2]; s[1] = xb[i2+1]; s[2] = xb[i2+2]; s[3] = xb[i2+3]; s[4] = xb[i2+4];
471:         PetscKernel_v_gets_A_times_w_5(xw,idiag,s);
472:         x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3]; x[i2+4] = xw[4];
473:         i2    -= 5;
474:         idiag -= 25;
475:         for (i=m-2; i>=0; i--) {
476:           v  = aa + 25*(diag[i]+1);
477:           vi = aj + diag[i] + 1;
478:           nz = ai[i+1] - diag[i] - 1;
479:           s[0] = xb[i2]; s[1] = xb[i2+1]; s[2] = xb[i2+2]; s[3] = xb[i2+3]; s[4] = xb[i2+4];
480:           while (nz--) {
481:             idx = 5*(*vi++);
482:             xw[0]  = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx]; xw[3] = x[3+idx]; xw[4] = x[4+idx];
483:             PetscKernel_v_gets_v_minus_A_times_w_5(s,v,xw);
484:             v  += 25;
485:           }
486:           PetscKernel_v_gets_A_times_w_5(xw,idiag,s);
487:           x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3]; x[i2+4] = xw[4];
488:           idiag  -= 25;
489:           i2     -= 5;
490:         }
491:         break;
492:       case 6:
493:         s[0]  = xb[i2]; s[1] = xb[i2+1]; s[2] = xb[i2+2]; s[3] = xb[i2+3]; s[4] = xb[i2+4]; s[5] = xb[i2+5];
494:         PetscKernel_v_gets_A_times_w_6(xw,idiag,s);
495:         x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3]; x[i2+4] = xw[4]; x[i2+5] = xw[5];
496:         i2    -= 6;
497:         idiag -= 36;
498:         for (i=m-2; i>=0; i--) {
499:           v  = aa + 36*(diag[i]+1);
500:           vi = aj + diag[i] + 1;
501:           nz = ai[i+1] - diag[i] - 1;
502:           s[0] = xb[i2]; s[1] = xb[i2+1]; s[2] = xb[i2+2]; s[3] = xb[i2+3]; s[4] = xb[i2+4]; s[5] = xb[i2+5];
503:           while (nz--) {
504:             idx = 6*(*vi++);
505:             xw[0] = x[idx];   xw[1] = x[1+idx]; xw[2] = x[2+idx];
506:             xw[3] = x[3+idx]; xw[4] = x[4+idx]; xw[5] = x[5+idx];
507:             PetscKernel_v_gets_v_minus_A_times_w_6(s,v,xw);
508:             v  += 36;
509:           }
510:           PetscKernel_v_gets_A_times_w_6(xw,idiag,s);
511:           x[i2] = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2]; x[i2+3] = xw[3]; x[i2+4] = xw[4]; x[i2+5] = xw[5];
512:           idiag  -= 36;
513:           i2     -= 6;
514:         }
515:         break;
516:       case 7:
517:         s[0] = xb[i2];   s[1] = xb[i2+1]; s[2] = xb[i2+2];
518:         s[3] = xb[i2+3]; s[4] = xb[i2+4]; s[5] = xb[i2+5]; s[6] = xb[i2+6];
519:         PetscKernel_v_gets_A_times_w_7(x,idiag,b);
520:         x[i2]   = xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2];
521:         x[i2+3] = xw[3]; x[i2+4] = xw[4]; x[i2+5] = xw[5]; x[i2+6] = xw[6];
522:         i2    -= 7;
523:         idiag -= 49;
524:         for (i=m-2; i>=0; i--) {
525:           v  = aa + 49*(diag[i]+1);
526:           vi = aj + diag[i] + 1;
527:           nz = ai[i+1] - diag[i] - 1;
528:           s[0] = xb[i2];   s[1] = xb[i2+1]; s[2] = xb[i2+2];
529:           s[3] = xb[i2+3]; s[4] = xb[i2+4]; s[5] = xb[i2+5]; s[6] = xb[i2+6];
530:           while (nz--) {
531:             idx = 7*(*vi++);
532:             xw[0] = x[idx];   xw[1] = x[1+idx]; xw[2] = x[2+idx];
533:             xw[3] = x[3+idx]; xw[4] = x[4+idx]; xw[5] = x[5+idx]; xw[6] = x[6+idx];
534:             PetscKernel_v_gets_v_minus_A_times_w_7(s,v,xw);
535:             v  += 49;
536:           }
537:           PetscKernel_v_gets_A_times_w_7(xw,idiag,s);
538:           x[i2] =   xw[0]; x[i2+1] = xw[1]; x[i2+2] = xw[2];
539:           x[i2+3] = xw[3]; x[i2+4] = xw[4]; x[i2+5] = xw[5]; x[i2+6] = xw[6];
540:           idiag  -= 49;
541:           i2     -= 7;
542:         }
543:         break;
544:       default:
545:         PetscArraycpy(w,xb+i2,bs);
546:         PetscKernel_w_gets_Ar_times_v(bs,bs,w,idiag,x+i2);
547:         i2    -= bs;
548:         idiag -= bs2;
549:         for (i=m-2; i>=0; i--) {
550:           v  = aa + bs2*(diag[i]+1);
551:           vi = aj + diag[i] + 1;
552:           nz = ai[i+1] - diag[i] - 1;

554:           PetscArraycpy(w,xb+i2,bs);
555:           /* copy all rows of x that are needed into contiguous space */
556:           workt = work;
557:           for (j=0; j<nz; j++) {
558:             PetscArraycpy(workt,x + bs*(*vi++),bs);
559:             workt += bs;
560:           }
561:           PetscKernel_w_gets_w_minus_Ar_times_v(bs,bs*nz,w,v,work);
562:           PetscKernel_w_gets_Ar_times_v(bs,bs,w,idiag,x+i2);

564:           idiag -= bs2;
565:           i2    -= bs;
566:         }
567:         break;
568:       }
569:       PetscLogFlops(1.0*bs2*(a->nz));
570:     }
571:     its--;
572:   }
573:   while (its--) {
574:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
575:       idiag = a->idiag;
576:       i2 = 0;
577:       switch (bs) {
578:       case 1:
579:         for (i=0; i<m; i++) {
580:           v  = aa + ai[i];
581:           vi = aj + ai[i];
582:           nz = ai[i+1] - ai[i];
583:           s[0] = b[i2];
584:           for (j=0; j<nz; j++) {
585:             xw[0] = x[vi[j]];
586:             PetscKernel_v_gets_v_minus_A_times_w_1(s,(v+j),xw);
587:           }
588:           PetscKernel_v_gets_A_times_w_1(xw,idiag,s);
589:           x[i2] += xw[0];
590:           idiag += 1;
591:           i2    += 1;
592:         }
593:         break;
594:       case 2:
595:         for (i=0; i<m; i++) {
596:           v  = aa + 4*ai[i];
597:           vi = aj + ai[i];
598:           nz = ai[i+1] - ai[i];
599:           s[0] = b[i2]; s[1] = b[i2+1];
600:           for (j=0; j<nz; j++) {
601:             idx = 2*vi[j];
602:             it  = 4*j;
603:             xw[0] = x[idx]; xw[1] = x[1+idx];
604:             PetscKernel_v_gets_v_minus_A_times_w_2(s,(v+it),xw);
605:           }
606:           PetscKernel_v_gets_A_times_w_2(xw,idiag,s);
607:           x[i2]  += xw[0]; x[i2+1] += xw[1];
608:           idiag  += 4;
609:           i2     += 2;
610:         }
611:         break;
612:       case 3:
613:         for (i=0; i<m; i++) {
614:           v  = aa + 9*ai[i];
615:           vi = aj + ai[i];
616:           nz = ai[i+1] - ai[i];
617:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2];
618:           while (nz--) {
619:             idx = 3*(*vi++);
620:             xw[0] = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx];
621:             PetscKernel_v_gets_v_minus_A_times_w_3(s,v,xw);
622:             v  += 9;
623:           }
624:           PetscKernel_v_gets_A_times_w_3(xw,idiag,s);
625:           x[i2] += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2];
626:           idiag  += 9;
627:           i2     += 3;
628:         }
629:         break;
630:       case 4:
631:         for (i=0; i<m; i++) {
632:           v  = aa + 16*ai[i];
633:           vi = aj + ai[i];
634:           nz = ai[i+1] - ai[i];
635:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3];
636:           while (nz--) {
637:             idx = 4*(*vi++);
638:             xw[0]  = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx]; xw[3] = x[3+idx];
639:             PetscKernel_v_gets_v_minus_A_times_w_4(s,v,xw);
640:             v  += 16;
641:           }
642:           PetscKernel_v_gets_A_times_w_4(xw,idiag,s);
643:           x[i2] += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2]; x[i2+3] += xw[3];
644:           idiag  += 16;
645:           i2     += 4;
646:         }
647:         break;
648:       case 5:
649:         for (i=0; i<m; i++) {
650:           v  = aa + 25*ai[i];
651:           vi = aj + ai[i];
652:           nz = ai[i+1] - ai[i];
653:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3]; s[4] = b[i2+4];
654:           while (nz--) {
655:             idx = 5*(*vi++);
656:             xw[0]  = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx]; xw[3] = x[3+idx]; xw[4] = x[4+idx];
657:             PetscKernel_v_gets_v_minus_A_times_w_5(s,v,xw);
658:             v  += 25;
659:           }
660:           PetscKernel_v_gets_A_times_w_5(xw,idiag,s);
661:           x[i2] += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2]; x[i2+3] += xw[3]; x[i2+4] += xw[4];
662:           idiag  += 25;
663:           i2     += 5;
664:         }
665:         break;
666:       case 6:
667:         for (i=0; i<m; i++) {
668:           v  = aa + 36*ai[i];
669:           vi = aj + ai[i];
670:           nz = ai[i+1] - ai[i];
671:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3]; s[4] = b[i2+4]; s[5] = b[i2+5];
672:           while (nz--) {
673:             idx = 6*(*vi++);
674:             xw[0] = x[idx];   xw[1] = x[1+idx]; xw[2] = x[2+idx];
675:             xw[3] = x[3+idx]; xw[4] = x[4+idx]; xw[5] = x[5+idx];
676:             PetscKernel_v_gets_v_minus_A_times_w_6(s,v,xw);
677:             v  += 36;
678:           }
679:           PetscKernel_v_gets_A_times_w_6(xw,idiag,s);
680:           x[i2] += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2];
681:           x[i2+3] += xw[3]; x[i2+4] += xw[4]; x[i2+5] += xw[5];
682:           idiag  += 36;
683:           i2     += 6;
684:         }
685:         break;
686:       case 7:
687:         for (i=0; i<m; i++) {
688:           v  = aa + 49*ai[i];
689:           vi = aj + ai[i];
690:           nz = ai[i+1] - ai[i];
691:           s[0] = b[i2];   s[1] = b[i2+1]; s[2] = b[i2+2];
692:           s[3] = b[i2+3]; s[4] = b[i2+4]; s[5] = b[i2+5]; s[6] = b[i2+6];
693:           while (nz--) {
694:             idx = 7*(*vi++);
695:             xw[0] = x[idx];   xw[1] = x[1+idx]; xw[2] = x[2+idx];
696:             xw[3] = x[3+idx]; xw[4] = x[4+idx]; xw[5] = x[5+idx]; xw[6] = x[6+idx];
697:             PetscKernel_v_gets_v_minus_A_times_w_7(s,v,xw);
698:             v  += 49;
699:           }
700:           PetscKernel_v_gets_A_times_w_7(xw,idiag,s);
701:           x[i2]   += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2];
702:           x[i2+3] += xw[3]; x[i2+4] += xw[4]; x[i2+5] += xw[5]; x[i2+6] += xw[6];
703:           idiag  += 49;
704:           i2     += 7;
705:         }
706:         break;
707:       default:
708:         for (i=0; i<m; i++) {
709:           v  = aa + bs2*ai[i];
710:           vi = aj + ai[i];
711:           nz = ai[i+1] - ai[i];

713:           PetscArraycpy(w,b+i2,bs);
714:           /* copy all rows of x that are needed into contiguous space */
715:           workt = work;
716:           for (j=0; j<nz; j++) {
717:             PetscArraycpy(workt,x + bs*(*vi++),bs);
718:             workt += bs;
719:           }
720:           PetscKernel_w_gets_w_minus_Ar_times_v(bs,bs*nz,w,v,work);
721:           PetscKernel_w_gets_w_plus_Ar_times_v(bs,bs,w,idiag,x+i2);

723:           idiag += bs2;
724:           i2    += bs;
725:         }
726:         break;
727:       }
728:       PetscLogFlops(2.0*bs2*a->nz);
729:     }
730:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
731:       idiag = a->idiag+bs2*(a->mbs-1);
732:       i2 = bs * (m-1);
733:       switch (bs) {
734:       case 1:
735:         for (i=m-1; i>=0; i--) {
736:           v  = aa + ai[i];
737:           vi = aj + ai[i];
738:           nz = ai[i+1] - ai[i];
739:           s[0] = b[i2];
740:           for (j=0; j<nz; j++) {
741:             xw[0] = x[vi[j]];
742:             PetscKernel_v_gets_v_minus_A_times_w_1(s,(v+j),xw);
743:           }
744:           PetscKernel_v_gets_A_times_w_1(xw,idiag,s);
745:           x[i2] += xw[0];
746:           idiag -= 1;
747:           i2    -= 1;
748:         }
749:         break;
750:       case 2:
751:         for (i=m-1; i>=0; i--) {
752:           v  = aa + 4*ai[i];
753:           vi = aj + ai[i];
754:           nz = ai[i+1] - ai[i];
755:           s[0] = b[i2]; s[1] = b[i2+1];
756:           for (j=0; j<nz; j++) {
757:             idx = 2*vi[j];
758:             it  = 4*j;
759:             xw[0] = x[idx]; xw[1] = x[1+idx];
760:             PetscKernel_v_gets_v_minus_A_times_w_2(s,(v+it),xw);
761:           }
762:           PetscKernel_v_gets_A_times_w_2(xw,idiag,s);
763:           x[i2]  += xw[0]; x[i2+1] += xw[1];
764:           idiag  -= 4;
765:           i2     -= 2;
766:         }
767:         break;
768:       case 3:
769:         for (i=m-1; i>=0; i--) {
770:           v  = aa + 9*ai[i];
771:           vi = aj + ai[i];
772:           nz = ai[i+1] - ai[i];
773:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2];
774:           while (nz--) {
775:             idx = 3*(*vi++);
776:             xw[0] = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx];
777:             PetscKernel_v_gets_v_minus_A_times_w_3(s,v,xw);
778:             v  += 9;
779:           }
780:           PetscKernel_v_gets_A_times_w_3(xw,idiag,s);
781:           x[i2] += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2];
782:           idiag  -= 9;
783:           i2     -= 3;
784:         }
785:         break;
786:       case 4:
787:         for (i=m-1; i>=0; i--) {
788:           v  = aa + 16*ai[i];
789:           vi = aj + ai[i];
790:           nz = ai[i+1] - ai[i];
791:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3];
792:           while (nz--) {
793:             idx = 4*(*vi++);
794:             xw[0]  = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx]; xw[3] = x[3+idx];
795:             PetscKernel_v_gets_v_minus_A_times_w_4(s,v,xw);
796:             v  += 16;
797:           }
798:           PetscKernel_v_gets_A_times_w_4(xw,idiag,s);
799:           x[i2] += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2]; x[i2+3] += xw[3];
800:           idiag  -= 16;
801:           i2     -= 4;
802:         }
803:         break;
804:       case 5:
805:         for (i=m-1; i>=0; i--) {
806:           v  = aa + 25*ai[i];
807:           vi = aj + ai[i];
808:           nz = ai[i+1] - ai[i];
809:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3]; s[4] = b[i2+4];
810:           while (nz--) {
811:             idx = 5*(*vi++);
812:             xw[0]  = x[idx]; xw[1] = x[1+idx]; xw[2] = x[2+idx]; xw[3] = x[3+idx]; xw[4] = x[4+idx];
813:             PetscKernel_v_gets_v_minus_A_times_w_5(s,v,xw);
814:             v  += 25;
815:           }
816:           PetscKernel_v_gets_A_times_w_5(xw,idiag,s);
817:           x[i2] += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2]; x[i2+3] += xw[3]; x[i2+4] += xw[4];
818:           idiag  -= 25;
819:           i2     -= 5;
820:         }
821:         break;
822:       case 6:
823:         for (i=m-1; i>=0; i--) {
824:           v  = aa + 36*ai[i];
825:           vi = aj + ai[i];
826:           nz = ai[i+1] - ai[i];
827:           s[0] = b[i2]; s[1] = b[i2+1]; s[2] = b[i2+2]; s[3] = b[i2+3]; s[4] = b[i2+4]; s[5] = b[i2+5];
828:           while (nz--) {
829:             idx = 6*(*vi++);
830:             xw[0] = x[idx];   xw[1] = x[1+idx]; xw[2] = x[2+idx];
831:             xw[3] = x[3+idx]; xw[4] = x[4+idx]; xw[5] = x[5+idx];
832:             PetscKernel_v_gets_v_minus_A_times_w_6(s,v,xw);
833:             v  += 36;
834:           }
835:           PetscKernel_v_gets_A_times_w_6(xw,idiag,s);
836:           x[i2] += xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2];
837:           x[i2+3] += xw[3]; x[i2+4] += xw[4]; x[i2+5] += xw[5];
838:           idiag  -= 36;
839:           i2     -= 6;
840:         }
841:         break;
842:       case 7:
843:         for (i=m-1; i>=0; i--) {
844:           v  = aa + 49*ai[i];
845:           vi = aj + ai[i];
846:           nz = ai[i+1] - ai[i];
847:           s[0] = b[i2];   s[1] = b[i2+1]; s[2] = b[i2+2];
848:           s[3] = b[i2+3]; s[4] = b[i2+4]; s[5] = b[i2+5]; s[6] = b[i2+6];
849:           while (nz--) {
850:             idx = 7*(*vi++);
851:             xw[0] = x[idx];   xw[1] = x[1+idx]; xw[2] = x[2+idx];
852:             xw[3] = x[3+idx]; xw[4] = x[4+idx]; xw[5] = x[5+idx]; xw[6] = x[6+idx];
853:             PetscKernel_v_gets_v_minus_A_times_w_7(s,v,xw);
854:             v  += 49;
855:           }
856:           PetscKernel_v_gets_A_times_w_7(xw,idiag,s);
857:           x[i2] +=   xw[0]; x[i2+1] += xw[1]; x[i2+2] += xw[2];
858:           x[i2+3] += xw[3]; x[i2+4] += xw[4]; x[i2+5] += xw[5]; x[i2+6] += xw[6];
859:           idiag  -= 49;
860:           i2     -= 7;
861:         }
862:         break;
863:       default:
864:         for (i=m-1; i>=0; i--) {
865:           v  = aa + bs2*ai[i];
866:           vi = aj + ai[i];
867:           nz = ai[i+1] - ai[i];

869:           PetscArraycpy(w,b+i2,bs);
870:           /* copy all rows of x that are needed into contiguous space */
871:           workt = work;
872:           for (j=0; j<nz; j++) {
873:             PetscArraycpy(workt,x + bs*(*vi++),bs);
874:             workt += bs;
875:           }
876:           PetscKernel_w_gets_w_minus_Ar_times_v(bs,bs*nz,w,v,work);
877:           PetscKernel_w_gets_w_plus_Ar_times_v(bs,bs,w,idiag,x+i2);

879:           idiag -= bs2;
880:           i2    -= bs;
881:         }
882:         break;
883:       }
884:       PetscLogFlops(2.0*bs2*(a->nz));
885:     }
886:   }
887:   VecRestoreArray(xx,&x);
888:   VecRestoreArrayRead(bb,&b);
889:   return(0);
890: }


893: /*
894:     Special version for direct calls from Fortran (Used in PETSc-fun3d)
895: */
896: #if defined(PETSC_HAVE_FORTRAN_CAPS)
897: #define matsetvaluesblocked4_ MATSETVALUESBLOCKED4
898: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
899: #define matsetvaluesblocked4_ matsetvaluesblocked4
900: #endif

902: PETSC_EXTERN void matsetvaluesblocked4_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[])
903: {
904:   Mat               A  = *AA;
905:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
906:   PetscInt          *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,N,m = *mm,n = *nn;
907:   PetscInt          *ai    =a->i,*ailen=a->ilen;
908:   PetscInt          *aj    =a->j,stepval,lastcol = -1;
909:   const PetscScalar *value = v;
910:   MatScalar         *ap,*aa = a->a,*bap;
911:   PetscErrorCode    ierr;

914:   if (A->rmap->bs != 4) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Can only be called with a block size of 4");
915:   stepval = (n-1)*4;
916:   for (k=0; k<m; k++) { /* loop over added rows */
917:     row  = im[k];
918:     rp   = aj + ai[row];
919:     ap   = aa + 16*ai[row];
920:     nrow = ailen[row];
921:     low  = 0;
922:     high = nrow;
923:     for (l=0; l<n; l++) { /* loop over added columns */
924:       col = in[l];
925:       if (col <= lastcol)  low = 0;
926:       else                high = nrow;
927:       lastcol = col;
928:       value   = v + k*(stepval+4 + l)*4;
929:       while (high-low > 7) {
930:         t = (low+high)/2;
931:         if (rp[t] > col) high = t;
932:         else             low  = t;
933:       }
934:       for (i=low; i<high; i++) {
935:         if (rp[i] > col) break;
936:         if (rp[i] == col) {
937:           bap = ap +  16*i;
938:           for (ii=0; ii<4; ii++,value+=stepval) {
939:             for (jj=ii; jj<16; jj+=4) {
940:               bap[jj] += *value++;
941:             }
942:           }
943:           goto noinsert2;
944:         }
945:       }
946:       N = nrow++ - 1;
947:       high++; /* added new column index thus must search to one higher than before */
948:       /* shift up all the later entries in this row */
949:       for (ii=N; ii>=i; ii--) {
950:         rp[ii+1] = rp[ii];
951:         PetscArraycpy(ap+16*(ii+1),ap+16*(ii),16);CHKERRV(ierr);
952:       }
953:       if (N >= i) {
954:         PetscArrayzero(ap+16*i,16);CHKERRV(ierr);
955:       }
956:       rp[i] = col;
957:       bap   = ap +  16*i;
958:       for (ii=0; ii<4; ii++,value+=stepval) {
959:         for (jj=ii; jj<16; jj+=4) {
960:           bap[jj] = *value++;
961:         }
962:       }
963:       noinsert2:;
964:       low = i;
965:     }
966:     ailen[row] = nrow;
967:   }
968:   PetscFunctionReturnVoid();
969: }

971: #if defined(PETSC_HAVE_FORTRAN_CAPS)
972: #define matsetvalues4_ MATSETVALUES4
973: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
974: #define matsetvalues4_ matsetvalues4
975: #endif

977: PETSC_EXTERN void matsetvalues4_(Mat *AA,PetscInt *mm,PetscInt *im,PetscInt *nn,PetscInt *in,PetscScalar *v)
978: {
979:   Mat         A  = *AA;
980:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
981:   PetscInt    *rp,k,low,high,t,row,nrow,i,col,l,N,n = *nn,m = *mm;
982:   PetscInt    *ai=a->i,*ailen=a->ilen;
983:   PetscInt    *aj=a->j,brow,bcol;
984:   PetscInt    ridx,cidx,lastcol = -1;
985:   MatScalar   *ap,value,*aa=a->a,*bap;

989:   for (k=0; k<m; k++) { /* loop over added rows */
990:     row  = im[k]; brow = row/4;
991:     rp   = aj + ai[brow];
992:     ap   = aa + 16*ai[brow];
993:     nrow = ailen[brow];
994:     low  = 0;
995:     high = nrow;
996:     for (l=0; l<n; l++) { /* loop over added columns */
997:       col   = in[l]; bcol = col/4;
998:       ridx  = row % 4; cidx = col % 4;
999:       value = v[l + k*n];
1000:       if (col <= lastcol)  low = 0;
1001:       else                high = nrow;
1002:       lastcol = col;
1003:       while (high-low > 7) {
1004:         t = (low+high)/2;
1005:         if (rp[t] > bcol) high = t;
1006:         else              low  = t;
1007:       }
1008:       for (i=low; i<high; i++) {
1009:         if (rp[i] > bcol) break;
1010:         if (rp[i] == bcol) {
1011:           bap   = ap +  16*i + 4*cidx + ridx;
1012:           *bap += value;
1013:           goto noinsert1;
1014:         }
1015:       }
1016:       N = nrow++ - 1;
1017:       high++; /* added new column thus must search to one higher than before */
1018:       /* shift up all the later entries in this row */
1019:       PetscArraymove(rp+i+1,rp+i,N-i+1);CHKERRV(ierr);
1020:       PetscArraymove(ap+16*i+16,ap+16*i,16*(N-i+1));CHKERRV(ierr);
1021:       PetscArrayzero(ap+16*i,16);CHKERRV(ierr);
1022:       rp[i]                    = bcol;
1023:       ap[16*i + 4*cidx + ridx] = value;
1024: noinsert1:;
1025:       low = i;
1026:     }
1027:     ailen[brow] = nrow;
1028:   }
1029:   PetscFunctionReturnVoid();
1030: }

1032: /*
1033:      Checks for missing diagonals
1034: */
1035: PetscErrorCode MatMissingDiagonal_SeqBAIJ(Mat A,PetscBool  *missing,PetscInt *d)
1036: {
1037:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
1039:   PetscInt       *diag,*ii = a->i,i;

1042:   MatMarkDiagonal_SeqBAIJ(A);
1043:   *missing = PETSC_FALSE;
1044:   if (A->rmap->n > 0 && !ii) {
1045:     *missing = PETSC_TRUE;
1046:     if (d) *d = 0;
1047:     PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");
1048:   } else {
1049:     PetscInt n;
1050:     n = PetscMin(a->mbs, a->nbs);
1051:     diag = a->diag;
1052:     for (i=0; i<n; i++) {
1053:       if (diag[i] >= ii[i+1]) {
1054:         *missing = PETSC_TRUE;
1055:         if (d) *d = i;
1056:         PetscInfo1(A,"Matrix is missing block diagonal number %D\n",i);
1057:         break;
1058:       }
1059:     }
1060:   }
1061:   return(0);
1062: }

1064: PetscErrorCode MatMarkDiagonal_SeqBAIJ(Mat A)
1065: {
1066:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
1068:   PetscInt       i,j,m = a->mbs;

1071:   if (!a->diag) {
1072:     PetscMalloc1(m,&a->diag);
1073:     PetscLogObjectMemory((PetscObject)A,m*sizeof(PetscInt));
1074:     a->free_diag = PETSC_TRUE;
1075:   }
1076:   for (i=0; i<m; i++) {
1077:     a->diag[i] = a->i[i+1];
1078:     for (j=a->i[i]; j<a->i[i+1]; j++) {
1079:       if (a->j[j] == i) {
1080:         a->diag[i] = j;
1081:         break;
1082:       }
1083:     }
1084:   }
1085:   return(0);
1086: }


1089: static PetscErrorCode MatGetRowIJ_SeqBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *inia[],const PetscInt *inja[],PetscBool  *done)
1090: {
1091:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
1093:   PetscInt       i,j,n = a->mbs,nz = a->i[n],*tia,*tja,bs = A->rmap->bs,k,l,cnt;
1094:   PetscInt       **ia = (PetscInt**)inia,**ja = (PetscInt**)inja;

1097:   *nn = n;
1098:   if (!ia) return(0);
1099:   if (symmetric) {
1100:     MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,PETSC_TRUE,0,0,&tia,&tja);
1101:     nz   = tia[n];
1102:   } else {
1103:     tia = a->i; tja = a->j;
1104:   }

1106:   if (!blockcompressed && bs > 1) {
1107:     (*nn) *= bs;
1108:     /* malloc & create the natural set of indices */
1109:     PetscMalloc1((n+1)*bs,ia);
1110:     if (n) {
1111:       (*ia)[0] = oshift;
1112:       for (j=1; j<bs; j++) {
1113:         (*ia)[j] = (tia[1]-tia[0])*bs+(*ia)[j-1];
1114:       }
1115:     }

1117:     for (i=1; i<n; i++) {
1118:       (*ia)[i*bs] = (tia[i]-tia[i-1])*bs + (*ia)[i*bs-1];
1119:       for (j=1; j<bs; j++) {
1120:         (*ia)[i*bs+j] = (tia[i+1]-tia[i])*bs + (*ia)[i*bs+j-1];
1121:       }
1122:     }
1123:     if (n) {
1124:       (*ia)[n*bs] = (tia[n]-tia[n-1])*bs + (*ia)[n*bs-1];
1125:     }

1127:     if (inja) {
1128:       PetscMalloc1(nz*bs*bs,ja);
1129:       cnt = 0;
1130:       for (i=0; i<n; i++) {
1131:         for (j=0; j<bs; j++) {
1132:           for (k=tia[i]; k<tia[i+1]; k++) {
1133:             for (l=0; l<bs; l++) {
1134:               (*ja)[cnt++] = bs*tja[k] + l;
1135:             }
1136:           }
1137:         }
1138:       }
1139:     }

1141:     if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
1142:       PetscFree(tia);
1143:       PetscFree(tja);
1144:     }
1145:   } else if (oshift == 1) {
1146:     if (symmetric) {
1147:       nz = tia[A->rmap->n/bs];
1148:       /*  add 1 to i and j indices */
1149:       for (i=0; i<A->rmap->n/bs+1; i++) tia[i] = tia[i] + 1;
1150:       *ia = tia;
1151:       if (ja) {
1152:         for (i=0; i<nz; i++) tja[i] = tja[i] + 1;
1153:         *ja = tja;
1154:       }
1155:     } else {
1156:       nz = a->i[A->rmap->n/bs];
1157:       /* malloc space and  add 1 to i and j indices */
1158:       PetscMalloc1(A->rmap->n/bs+1,ia);
1159:       for (i=0; i<A->rmap->n/bs+1; i++) (*ia)[i] = a->i[i] + 1;
1160:       if (ja) {
1161:         PetscMalloc1(nz,ja);
1162:         for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1;
1163:       }
1164:     }
1165:   } else {
1166:     *ia = tia;
1167:     if (ja) *ja = tja;
1168:   }
1169:   return(0);
1170: }

1172: static PetscErrorCode MatRestoreRowIJ_SeqBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool  *done)
1173: {

1177:   if (!ia) return(0);
1178:   if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
1179:     PetscFree(*ia);
1180:     if (ja) {PetscFree(*ja);}
1181:   }
1182:   return(0);
1183: }

1185: PetscErrorCode MatDestroy_SeqBAIJ(Mat A)
1186: {
1187:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;

1191: #if defined(PETSC_USE_LOG)
1192:   PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->N,A->cmap->n,a->nz);
1193: #endif
1194:   MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);
1195:   ISDestroy(&a->row);
1196:   ISDestroy(&a->col);
1197:   if (a->free_diag) {PetscFree(a->diag);}
1198:   PetscFree(a->idiag);
1199:   if (a->free_imax_ilen) {PetscFree2(a->imax,a->ilen);}
1200:   PetscFree(a->solve_work);
1201:   PetscFree(a->mult_work);
1202:   PetscFree(a->sor_workt);
1203:   PetscFree(a->sor_work);
1204:   ISDestroy(&a->icol);
1205:   PetscFree(a->saved_values);
1206:   PetscFree2(a->compressedrow.i,a->compressedrow.rindex);

1208:   MatDestroy(&a->sbaijMat);
1209:   MatDestroy(&a->parent);
1210:   PetscFree(A->data);

1212:   PetscObjectChangeTypeName((PetscObject)A,0);
1213:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBAIJGetArray_C",NULL);
1214:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBAIJRestoreArray_C",NULL);
1215:   PetscObjectComposeFunction((PetscObject)A,"MatInvertBlockDiagonal_C",NULL);
1216:   PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);
1217:   PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);
1218:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBAIJSetColumnIndices_C",NULL);
1219:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_seqaij_C",NULL);
1220:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_seqsbaij_C",NULL);
1221:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",NULL);
1222:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBAIJSetPreallocationCSR_C",NULL);
1223:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_seqbstrm_C",NULL);
1224:   PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);
1225: #if defined(PETSC_HAVE_HYPRE)
1226:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_hypre_C",NULL);
1227: #endif
1228:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_is_C",NULL);
1229:   return(0);
1230: }

1232: PetscErrorCode MatSetOption_SeqBAIJ(Mat A,MatOption op,PetscBool flg)
1233: {
1234:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;

1238:   switch (op) {
1239:   case MAT_ROW_ORIENTED:
1240:     a->roworiented = flg;
1241:     break;
1242:   case MAT_KEEP_NONZERO_PATTERN:
1243:     a->keepnonzeropattern = flg;
1244:     break;
1245:   case MAT_NEW_NONZERO_LOCATIONS:
1246:     a->nonew = (flg ? 0 : 1);
1247:     break;
1248:   case MAT_NEW_NONZERO_LOCATION_ERR:
1249:     a->nonew = (flg ? -1 : 0);
1250:     break;
1251:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1252:     a->nonew = (flg ? -2 : 0);
1253:     break;
1254:   case MAT_UNUSED_NONZERO_LOCATION_ERR:
1255:     a->nounused = (flg ? -1 : 0);
1256:     break;
1257:   case MAT_NEW_DIAGONALS:
1258:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1259:   case MAT_USE_HASH_TABLE:
1260:   case MAT_SORTED_FULL:
1261:     PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);
1262:     break;
1263:   case MAT_SPD:
1264:   case MAT_SYMMETRIC:
1265:   case MAT_STRUCTURALLY_SYMMETRIC:
1266:   case MAT_HERMITIAN:
1267:   case MAT_SYMMETRY_ETERNAL:
1268:   case MAT_SUBMAT_SINGLEIS:
1269:   case MAT_STRUCTURE_ONLY:
1270:     /* These options are handled directly by MatSetOption() */
1271:     break;
1272:   default:
1273:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op);
1274:   }
1275:   return(0);
1276: }

1278: /* used for both SeqBAIJ and SeqSBAIJ matrices */
1279: PetscErrorCode MatGetRow_SeqBAIJ_private(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v,PetscInt *ai,PetscInt *aj,PetscScalar *aa)
1280: {
1282:   PetscInt       itmp,i,j,k,M,bn,bp,*idx_i,bs,bs2;
1283:   MatScalar      *aa_i;
1284:   PetscScalar    *v_i;

1287:   bs  = A->rmap->bs;
1288:   bs2 = bs*bs;
1289:   if (row < 0 || row >= A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range", row);

1291:   bn  = row/bs;   /* Block number */
1292:   bp  = row % bs; /* Block Position */
1293:   M   = ai[bn+1] - ai[bn];
1294:   *nz = bs*M;

1296:   if (v) {
1297:     *v = 0;
1298:     if (*nz) {
1299:       PetscMalloc1(*nz,v);
1300:       for (i=0; i<M; i++) { /* for each block in the block row */
1301:         v_i  = *v + i*bs;
1302:         aa_i = aa + bs2*(ai[bn] + i);
1303:         for (j=bp,k=0; j<bs2; j+=bs,k++) v_i[k] = aa_i[j];
1304:       }
1305:     }
1306:   }

1308:   if (idx) {
1309:     *idx = 0;
1310:     if (*nz) {
1311:       PetscMalloc1(*nz,idx);
1312:       for (i=0; i<M; i++) { /* for each block in the block row */
1313:         idx_i = *idx + i*bs;
1314:         itmp  = bs*aj[ai[bn] + i];
1315:         for (j=0; j<bs; j++) idx_i[j] = itmp++;
1316:       }
1317:     }
1318:   }
1319:   return(0);
1320: }

1322: PetscErrorCode MatGetRow_SeqBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1323: {
1324:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;

1328:   MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a);
1329:   return(0);
1330: }

1332: PetscErrorCode MatRestoreRow_SeqBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1333: {

1337:   if (idx) {PetscFree(*idx);}
1338:   if (v)   {PetscFree(*v);}
1339:   return(0);
1340: }

1342: PetscErrorCode MatTranspose_SeqBAIJ(Mat A,MatReuse reuse,Mat *B)
1343: {
1344:   Mat_SeqBAIJ    *a=(Mat_SeqBAIJ*)A->data,*at;
1345:   Mat            C;
1347:   PetscInt       i,j,k,*aj=a->j,*ai=a->i,bs=A->rmap->bs,mbs=a->mbs,nbs=a->nbs,*atfill;
1348:   PetscInt       bs2=a->bs2,*ati,*atj,anzj,kr;
1349:   MatScalar      *ata,*aa=a->a;

1352:   PetscCalloc1(1+nbs,&atfill);
1353:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1354:     for (i=0; i<ai[mbs]; i++) atfill[aj[i]] += 1; /* count num of non-zeros in row aj[i] */

1356:     MatCreate(PetscObjectComm((PetscObject)A),&C);
1357:     MatSetSizes(C,A->cmap->n,A->rmap->N,A->cmap->n,A->rmap->N);
1358:     MatSetType(C,((PetscObject)A)->type_name);
1359:     MatSeqBAIJSetPreallocation(C,bs,0,atfill);

1361:     at  = (Mat_SeqBAIJ*)C->data;
1362:     ati = at->i;
1363:     for (i=0; i<nbs; i++) at->ilen[i] = at->imax[i] = ati[i+1] - ati[i];
1364:   } else {
1365:     C = *B;
1366:     at = (Mat_SeqBAIJ*)C->data;
1367:     ati = at->i;
1368:   }

1370:   atj = at->j;
1371:   ata = at->a;

1373:   /* Copy ati into atfill so we have locations of the next free space in atj */
1374:   PetscArraycpy(atfill,ati,nbs);

1376:   /* Walk through A row-wise and mark nonzero entries of A^T. */
1377:   for (i=0; i<mbs; i++) {
1378:     anzj = ai[i+1] - ai[i];
1379:     for (j=0; j<anzj; j++) {
1380:       atj[atfill[*aj]] = i;
1381:       for (kr=0; kr<bs; kr++) {
1382:         for (k=0; k<bs; k++) {
1383:           ata[bs2*atfill[*aj]+k*bs+kr] = *aa++;
1384:         }
1385:       }
1386:       atfill[*aj++] += 1;
1387:     }
1388:   }
1389:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1390:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

1392:   /* Clean up temporary space and complete requests. */
1393:   PetscFree(atfill);

1395:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1396:     MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
1397:     *B = C;
1398:   } else {
1399:     MatHeaderMerge(A,&C);
1400:   }
1401:   return(0);
1402: }

1404: PetscErrorCode MatIsTranspose_SeqBAIJ(Mat A,Mat B,PetscReal tol,PetscBool  *f)
1405: {
1407:   Mat            Btrans;

1410:   *f   = PETSC_FALSE;
1411:   MatTranspose_SeqBAIJ(A,MAT_INITIAL_MATRIX,&Btrans);
1412:   MatEqual_SeqBAIJ(B,Btrans,f);
1413:   MatDestroy(&Btrans);
1414:   return(0);
1415: }

1417: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
1418: PetscErrorCode MatView_SeqBAIJ_Binary(Mat mat,PetscViewer viewer)
1419: {
1420:   Mat_SeqBAIJ    *A = (Mat_SeqBAIJ*)mat->data;
1421:   PetscInt       header[4],M,N,m,bs,nz,cnt,i,j,k,l;
1422:   PetscInt       *rowlens,*colidxs;
1423:   PetscScalar    *matvals;

1427:   PetscViewerSetUp(viewer);

1429:   M  = mat->rmap->N;
1430:   N  = mat->cmap->N;
1431:   m  = mat->rmap->n;
1432:   bs = mat->rmap->bs;
1433:   nz = bs*bs*A->nz;

1435:   /* write matrix header */
1436:   header[0] = MAT_FILE_CLASSID;
1437:   header[1] = M; header[2] = N; header[3] = nz;
1438:   PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);

1440:   /* store row lengths */
1441:   PetscMalloc1(m,&rowlens);
1442:   for (cnt=0, i=0; i<A->mbs; i++)
1443:     for (j=0; j<bs; j++)
1444:       rowlens[cnt++] = bs*(A->i[i+1] - A->i[i]);
1445:   PetscViewerBinaryWrite(viewer,rowlens,m,PETSC_INT);
1446:   PetscFree(rowlens);

1448:   /* store column indices  */
1449:   PetscMalloc1(nz,&colidxs);
1450:   for (cnt=0, i=0; i<A->mbs; i++)
1451:     for (k=0; k<bs; k++)
1452:       for (j=A->i[i]; j<A->i[i+1]; j++)
1453:         for (l=0; l<bs; l++)
1454:           colidxs[cnt++] = bs*A->j[j] + l;
1455:   if (cnt != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Internal PETSc error: cnt = %D nz = %D",cnt,nz);
1456:   PetscViewerBinaryWrite(viewer,colidxs,nz,PETSC_INT);
1457:   PetscFree(colidxs);

1459:   /* store nonzero values */
1460:   PetscMalloc1(nz,&matvals);
1461:   for (cnt=0, i=0; i<A->mbs; i++)
1462:     for (k=0; k<bs; k++)
1463:       for (j=A->i[i]; j<A->i[i+1]; j++)
1464:         for (l=0; l<bs; l++)
1465:           matvals[cnt++] = A->a[bs*(bs*j + l) + k];
1466:   if (cnt != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Internal PETSc error: cnt = %D nz = %D",cnt,nz);
1467:   PetscViewerBinaryWrite(viewer,matvals,nz,PETSC_SCALAR);
1468:   PetscFree(matvals);

1470:   /* write block size option to the viewer's .info file */
1471:   MatView_Binary_BlockSizes(mat,viewer);
1472:   return(0);
1473: }

1475: static PetscErrorCode MatView_SeqBAIJ_ASCII_structonly(Mat A,PetscViewer viewer)
1476: {
1478:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
1479:   PetscInt       i,bs = A->rmap->bs,k;

1482:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1483:   for (i=0; i<a->mbs; i++) {
1484:     PetscViewerASCIIPrintf(viewer,"row %D-%D:",i*bs,i*bs+bs-1);
1485:     for (k=a->i[i]; k<a->i[i+1]; k++) {
1486:       PetscViewerASCIIPrintf(viewer," (%D-%D) ",bs*a->j[k],bs*a->j[k]+bs-1);
1487:     }
1488:     PetscViewerASCIIPrintf(viewer,"\n");
1489:   }
1490:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1491:   return(0);
1492: }

1494: static PetscErrorCode MatView_SeqBAIJ_ASCII(Mat A,PetscViewer viewer)
1495: {
1496:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
1497:   PetscErrorCode    ierr;
1498:   PetscInt          i,j,bs = A->rmap->bs,k,l,bs2=a->bs2;
1499:   PetscViewerFormat format;

1502:   if (A->structure_only) {
1503:     MatView_SeqBAIJ_ASCII_structonly(A,viewer);
1504:     return(0);
1505:   }

1507:   PetscViewerGetFormat(viewer,&format);
1508:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1509:     PetscViewerASCIIPrintf(viewer,"  block size is %D\n",bs);
1510:   } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
1511:     const char *matname;
1512:     Mat        aij;
1513:     MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);
1514:     PetscObjectGetName((PetscObject)A,&matname);
1515:     PetscObjectSetName((PetscObject)aij,matname);
1516:     MatView(aij,viewer);
1517:     MatDestroy(&aij);
1518:   } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
1519:       return(0);
1520:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1521:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1522:     for (i=0; i<a->mbs; i++) {
1523:       for (j=0; j<bs; j++) {
1524:         PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
1525:         for (k=a->i[i]; k<a->i[i+1]; k++) {
1526:           for (l=0; l<bs; l++) {
1527: #if defined(PETSC_USE_COMPLEX)
1528:             if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
1529:               PetscViewerASCIIPrintf(viewer," (%D, %g + %gi) ",bs*a->j[k]+l,
1530:                                             (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
1531:             } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
1532:               PetscViewerASCIIPrintf(viewer," (%D, %g - %gi) ",bs*a->j[k]+l,
1533:                                             (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
1534:             } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) {
1535:               PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
1536:             }
1537: #else
1538:             if (a->a[bs2*k + l*bs + j] != 0.0) {
1539:               PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
1540:             }
1541: #endif
1542:           }
1543:         }
1544:         PetscViewerASCIIPrintf(viewer,"\n");
1545:       }
1546:     }
1547:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1548:   } else {
1549:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1550:     for (i=0; i<a->mbs; i++) {
1551:       for (j=0; j<bs; j++) {
1552:         PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);
1553:         for (k=a->i[i]; k<a->i[i+1]; k++) {
1554:           for (l=0; l<bs; l++) {
1555: #if defined(PETSC_USE_COMPLEX)
1556:             if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) {
1557:               PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l,
1558:                                             (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
1559:             } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) {
1560:               PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l,
1561:                                             (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));
1562:             } else {
1563:               PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));
1564:             }
1565: #else
1566:             PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);
1567: #endif
1568:           }
1569:         }
1570:         PetscViewerASCIIPrintf(viewer,"\n");
1571:       }
1572:     }
1573:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1574:   }
1575:   PetscViewerFlush(viewer);
1576:   return(0);
1577: }

1579:  #include <petscdraw.h>
1580: static PetscErrorCode MatView_SeqBAIJ_Draw_Zoom(PetscDraw draw,void *Aa)
1581: {
1582:   Mat               A = (Mat) Aa;
1583:   Mat_SeqBAIJ       *a=(Mat_SeqBAIJ*)A->data;
1584:   PetscErrorCode    ierr;
1585:   PetscInt          row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2;
1586:   PetscReal         xl,yl,xr,yr,x_l,x_r,y_l,y_r;
1587:   MatScalar         *aa;
1588:   PetscViewer       viewer;
1589:   PetscViewerFormat format;

1592:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1593:   PetscViewerGetFormat(viewer,&format);
1594:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

1596:   /* loop over matrix elements drawing boxes */

1598:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1599:     PetscDrawCollectiveBegin(draw);
1600:     /* Blue for negative, Cyan for zero and  Red for positive */
1601:     color = PETSC_DRAW_BLUE;
1602:     for (i=0,row=0; i<mbs; i++,row+=bs) {
1603:       for (j=a->i[i]; j<a->i[i+1]; j++) {
1604:         y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
1605:         x_l = a->j[j]*bs; x_r = x_l + 1.0;
1606:         aa  = a->a + j*bs2;
1607:         for (k=0; k<bs; k++) {
1608:           for (l=0; l<bs; l++) {
1609:             if (PetscRealPart(*aa++) >=  0.) continue;
1610:             PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
1611:           }
1612:         }
1613:       }
1614:     }
1615:     color = PETSC_DRAW_CYAN;
1616:     for (i=0,row=0; i<mbs; i++,row+=bs) {
1617:       for (j=a->i[i]; j<a->i[i+1]; j++) {
1618:         y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
1619:         x_l = a->j[j]*bs; x_r = x_l + 1.0;
1620:         aa  = a->a + j*bs2;
1621:         for (k=0; k<bs; k++) {
1622:           for (l=0; l<bs; l++) {
1623:             if (PetscRealPart(*aa++) != 0.) continue;
1624:             PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
1625:           }
1626:         }
1627:       }
1628:     }
1629:     color = PETSC_DRAW_RED;
1630:     for (i=0,row=0; i<mbs; i++,row+=bs) {
1631:       for (j=a->i[i]; j<a->i[i+1]; j++) {
1632:         y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
1633:         x_l = a->j[j]*bs; x_r = x_l + 1.0;
1634:         aa  = a->a + j*bs2;
1635:         for (k=0; k<bs; k++) {
1636:           for (l=0; l<bs; l++) {
1637:             if (PetscRealPart(*aa++) <= 0.) continue;
1638:             PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
1639:           }
1640:         }
1641:       }
1642:     }
1643:     PetscDrawCollectiveEnd(draw);
1644:   } else {
1645:     /* use contour shading to indicate magnitude of values */
1646:     /* first determine max of all nonzero values */
1647:     PetscReal minv = 0.0, maxv = 0.0;
1648:     PetscDraw popup;

1650:     for (i=0; i<a->nz*a->bs2; i++) {
1651:       if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]);
1652:     }
1653:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1654:     PetscDrawGetPopup(draw,&popup);
1655:     PetscDrawScalePopup(popup,0.0,maxv);

1657:     PetscDrawCollectiveBegin(draw);
1658:     for (i=0,row=0; i<mbs; i++,row+=bs) {
1659:       for (j=a->i[i]; j<a->i[i+1]; j++) {
1660:         y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0;
1661:         x_l = a->j[j]*bs; x_r = x_l + 1.0;
1662:         aa  = a->a + j*bs2;
1663:         for (k=0; k<bs; k++) {
1664:           for (l=0; l<bs; l++) {
1665:             MatScalar v = *aa++;
1666:             color = PetscDrawRealToColor(PetscAbsScalar(v),minv,maxv);
1667:             PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);
1668:           }
1669:         }
1670:       }
1671:     }
1672:     PetscDrawCollectiveEnd(draw);
1673:   }
1674:   return(0);
1675: }

1677: static PetscErrorCode MatView_SeqBAIJ_Draw(Mat A,PetscViewer viewer)
1678: {
1680:   PetscReal      xl,yl,xr,yr,w,h;
1681:   PetscDraw      draw;
1682:   PetscBool      isnull;

1685:   PetscViewerDrawGetDraw(viewer,0,&draw);
1686:   PetscDrawIsNull(draw,&isnull);
1687:   if (isnull) return(0);

1689:   xr   = A->cmap->n; yr = A->rmap->N; h = yr/10.0; w = xr/10.0;
1690:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1691:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1692:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1693:   PetscDrawZoom(draw,MatView_SeqBAIJ_Draw_Zoom,A);
1694:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1695:   PetscDrawSave(draw);
1696:   return(0);
1697: }

1699: PetscErrorCode MatView_SeqBAIJ(Mat A,PetscViewer viewer)
1700: {
1702:   PetscBool      iascii,isbinary,isdraw;

1705:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1706:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1707:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1708:   if (iascii) {
1709:     MatView_SeqBAIJ_ASCII(A,viewer);
1710:   } else if (isbinary) {
1711:     MatView_SeqBAIJ_Binary(A,viewer);
1712:   } else if (isdraw) {
1713:     MatView_SeqBAIJ_Draw(A,viewer);
1714:   } else {
1715:     Mat B;
1716:     MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
1717:     MatView(B,viewer);
1718:     MatDestroy(&B);
1719:   }
1720:   return(0);
1721: }


1724: PetscErrorCode MatGetValues_SeqBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
1725: {
1726:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
1727:   PetscInt    *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
1728:   PetscInt    *ai = a->i,*ailen = a->ilen;
1729:   PetscInt    brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2;
1730:   MatScalar   *ap,*aa = a->a;

1733:   for (k=0; k<m; k++) { /* loop over rows */
1734:     row = im[k]; brow = row/bs;
1735:     if (row < 0) {v += n; continue;} /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row"); */
1736:     if (row >= A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D too large", row);
1737:     rp   = aj + ai[brow]; ap = aa + bs2*ai[brow];
1738:     nrow = ailen[brow];
1739:     for (l=0; l<n; l++) { /* loop over columns */
1740:       if (in[l] < 0) {v++; continue;} /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column"); */
1741:       if (in[l] >= A->cmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D too large", in[l]);
1742:       col  = in[l];
1743:       bcol = col/bs;
1744:       cidx = col%bs;
1745:       ridx = row%bs;
1746:       high = nrow;
1747:       low  = 0; /* assume unsorted */
1748:       while (high-low > 5) {
1749:         t = (low+high)/2;
1750:         if (rp[t] > bcol) high = t;
1751:         else             low  = t;
1752:       }
1753:       for (i=low; i<high; i++) {
1754:         if (rp[i] > bcol) break;
1755:         if (rp[i] == bcol) {
1756:           *v++ = ap[bs2*i+bs*cidx+ridx];
1757:           goto finished;
1758:         }
1759:       }
1760:       *v++ = 0.0;
1761: finished:;
1762:     }
1763:   }
1764:   return(0);
1765: }

1767: PetscErrorCode MatSetValuesBlocked_SeqBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
1768: {
1769:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
1770:   PetscInt          *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1;
1771:   PetscInt          *imax=a->imax,*ai=a->i,*ailen=a->ilen;
1772:   PetscErrorCode    ierr;
1773:   PetscInt          *aj        =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval;
1774:   PetscBool         roworiented=a->roworiented;
1775:   const PetscScalar *value     = v;
1776:   MatScalar         *ap=NULL,*aa = a->a,*bap;

1779:   if (roworiented) {
1780:     stepval = (n-1)*bs;
1781:   } else {
1782:     stepval = (m-1)*bs;
1783:   }
1784:   for (k=0; k<m; k++) { /* loop over added rows */
1785:     row = im[k];
1786:     if (row < 0) continue;
1787: #if defined(PETSC_USE_DEBUG)
1788:     if (row >= a->mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block row index too large %D max %D",row,a->mbs-1);
1789: #endif
1790:     rp   = aj + ai[row];
1791:     if (!A->structure_only) ap = aa + bs2*ai[row];
1792:     rmax = imax[row];
1793:     nrow = ailen[row];
1794:     low  = 0;
1795:     high = nrow;
1796:     for (l=0; l<n; l++) { /* loop over added columns */
1797:       if (in[l] < 0) continue;
1798: #if defined(PETSC_USE_DEBUG)
1799:       if (in[l] >= a->nbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block column index too large %D max %D",in[l],a->nbs-1);
1800: #endif
1801:       col = in[l];
1802:       if (!A->structure_only) {
1803:         if (roworiented) {
1804:           value = v + (k*(stepval+bs) + l)*bs;
1805:         } else {
1806:           value = v + (l*(stepval+bs) + k)*bs;
1807:         }
1808:       }
1809:       if (col <= lastcol) low = 0;
1810:       else high = nrow;
1811:       lastcol = col;
1812:       while (high-low > 7) {
1813:         t = (low+high)/2;
1814:         if (rp[t] > col) high = t;
1815:         else             low  = t;
1816:       }
1817:       for (i=low; i<high; i++) {
1818:         if (rp[i] > col) break;
1819:         if (rp[i] == col) {
1820:           if (A->structure_only) goto noinsert2;
1821:           bap = ap +  bs2*i;
1822:           if (roworiented) {
1823:             if (is == ADD_VALUES) {
1824:               for (ii=0; ii<bs; ii++,value+=stepval) {
1825:                 for (jj=ii; jj<bs2; jj+=bs) {
1826:                   bap[jj] += *value++;
1827:                 }
1828:               }
1829:             } else {
1830:               for (ii=0; ii<bs; ii++,value+=stepval) {
1831:                 for (jj=ii; jj<bs2; jj+=bs) {
1832:                   bap[jj] = *value++;
1833:                 }
1834:               }
1835:             }
1836:           } else {
1837:             if (is == ADD_VALUES) {
1838:               for (ii=0; ii<bs; ii++,value+=bs+stepval) {
1839:                 for (jj=0; jj<bs; jj++) {
1840:                   bap[jj] += value[jj];
1841:                 }
1842:                 bap += bs;
1843:               }
1844:             } else {
1845:               for (ii=0; ii<bs; ii++,value+=bs+stepval) {
1846:                 for (jj=0; jj<bs; jj++) {
1847:                   bap[jj]  = value[jj];
1848:                 }
1849:                 bap += bs;
1850:               }
1851:             }
1852:           }
1853:           goto noinsert2;
1854:         }
1855:       }
1856:       if (nonew == 1) goto noinsert2;
1857:       if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new blocked index new nonzero block (%D, %D) in the matrix", row, col);
1858:       if (A->structure_only) {
1859:         MatSeqXAIJReallocateAIJ_structure_only(A,a->mbs,bs2,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar);
1860:       } else {
1861:         MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
1862:       }
1863:       N = nrow++ - 1; high++;
1864:       /* shift up all the later entries in this row */
1865:       PetscArraymove(rp+i+1,rp+i,N-i+1);
1866:       rp[i] = col;
1867:       if (!A->structure_only) {
1868:         PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
1869:         bap   = ap +  bs2*i;
1870:         if (roworiented) {
1871:           for (ii=0; ii<bs; ii++,value+=stepval) {
1872:             for (jj=ii; jj<bs2; jj+=bs) {
1873:               bap[jj] = *value++;
1874:             }
1875:           }
1876:         } else {
1877:           for (ii=0; ii<bs; ii++,value+=stepval) {
1878:             for (jj=0; jj<bs; jj++) {
1879:               *bap++ = *value++;
1880:             }
1881:           }
1882:         }
1883:       }
1884: noinsert2:;
1885:       low = i;
1886:     }
1887:     ailen[row] = nrow;
1888:   }
1889:   return(0);
1890: }

1892: PetscErrorCode MatAssemblyEnd_SeqBAIJ(Mat A,MatAssemblyType mode)
1893: {
1894:   Mat_SeqBAIJ    *a     = (Mat_SeqBAIJ*)A->data;
1895:   PetscInt       fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
1896:   PetscInt       m      = A->rmap->N,*ip,N,*ailen = a->ilen;
1898:   PetscInt       mbs  = a->mbs,bs2 = a->bs2,rmax = 0;
1899:   MatScalar      *aa  = a->a,*ap;
1900:   PetscReal      ratio=0.6;

1903:   if (mode == MAT_FLUSH_ASSEMBLY) return(0);

1905:   if (m) rmax = ailen[0];
1906:   for (i=1; i<mbs; i++) {
1907:     /* move each row back by the amount of empty slots (fshift) before it*/
1908:     fshift += imax[i-1] - ailen[i-1];
1909:     rmax    = PetscMax(rmax,ailen[i]);
1910:     if (fshift) {
1911:       ip = aj + ai[i];
1912:       ap = aa + bs2*ai[i];
1913:       N  = ailen[i];
1914:       PetscArraymove(ip-fshift,ip,N);
1915:       if (!A->structure_only) {
1916:         PetscArraymove(ap-bs2*fshift,ap,bs2*N);
1917:       }
1918:     }
1919:     ai[i] = ai[i-1] + ailen[i-1];
1920:   }
1921:   if (mbs) {
1922:     fshift += imax[mbs-1] - ailen[mbs-1];
1923:     ai[mbs] = ai[mbs-1] + ailen[mbs-1];
1924:   }

1926:   /* reset ilen and imax for each row */
1927:   a->nonzerorowcnt = 0;
1928:   if (A->structure_only) {
1929:     PetscFree2(a->imax,a->ilen);
1930:   } else { /* !A->structure_only */
1931:     for (i=0; i<mbs; i++) {
1932:       ailen[i] = imax[i] = ai[i+1] - ai[i];
1933:       a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0);
1934:     }
1935:   }
1936:   a->nz = ai[mbs];

1938:   /* diagonals may have moved, so kill the diagonal pointers */
1939:   a->idiagvalid = PETSC_FALSE;
1940:   if (fshift && a->diag) {
1941:     PetscFree(a->diag);
1942:     PetscLogObjectMemory((PetscObject)A,-(mbs+1)*sizeof(PetscInt));
1943:     a->diag = 0;
1944:   }
1945:   if (fshift && a->nounused == -1) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D block size %D, %D unneeded", m, A->cmap->n, A->rmap->bs, fshift*bs2);
1946:   PetscInfo5(A,"Matrix size: %D X %D, block size %D; storage space: %D unneeded, %D used\n",m,A->cmap->n,A->rmap->bs,fshift*bs2,a->nz*bs2);
1947:   PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);
1948:   PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);

1950:   A->info.mallocs    += a->reallocs;
1951:   a->reallocs         = 0;
1952:   A->info.nz_unneeded = (PetscReal)fshift*bs2;
1953:   a->rmax             = rmax;

1955:   if (!A->structure_only) {
1956:     MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,mbs,ratio);
1957:   }
1958:   return(0);
1959: }

1961: /*
1962:    This function returns an array of flags which indicate the locations of contiguous
1963:    blocks that should be zeroed. for eg: if bs = 3  and is = [0,1,2,3,5,6,7,8,9]
1964:    then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)]
1965:    Assume: sizes should be long enough to hold all the values.
1966: */
1967: static PetscErrorCode MatZeroRows_SeqBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max)
1968: {
1969:   PetscInt  i,j,k,row;
1970:   PetscBool flg;

1973:   for (i=0,j=0; i<n; j++) {
1974:     row = idx[i];
1975:     if (row%bs!=0) { /* Not the begining of a block */
1976:       sizes[j] = 1;
1977:       i++;
1978:     } else if (i+bs > n) { /* complete block doesn't exist (at idx end) */
1979:       sizes[j] = 1;         /* Also makes sure atleast 'bs' values exist for next else */
1980:       i++;
1981:     } else { /* Begining of the block, so check if the complete block exists */
1982:       flg = PETSC_TRUE;
1983:       for (k=1; k<bs; k++) {
1984:         if (row+k != idx[i+k]) { /* break in the block */
1985:           flg = PETSC_FALSE;
1986:           break;
1987:         }
1988:       }
1989:       if (flg) { /* No break in the bs */
1990:         sizes[j] = bs;
1991:         i       += bs;
1992:       } else {
1993:         sizes[j] = 1;
1994:         i++;
1995:       }
1996:     }
1997:   }
1998:   *bs_max = j;
1999:   return(0);
2000: }

2002: PetscErrorCode MatZeroRows_SeqBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
2003: {
2004:   Mat_SeqBAIJ       *baij=(Mat_SeqBAIJ*)A->data;
2005:   PetscErrorCode    ierr;
2006:   PetscInt          i,j,k,count,*rows;
2007:   PetscInt          bs=A->rmap->bs,bs2=baij->bs2,*sizes,row,bs_max;
2008:   PetscScalar       zero = 0.0;
2009:   MatScalar         *aa;
2010:   const PetscScalar *xx;
2011:   PetscScalar       *bb;

2014:   /* fix right hand side if needed */
2015:   if (x && b) {
2016:     VecGetArrayRead(x,&xx);
2017:     VecGetArray(b,&bb);
2018:     for (i=0; i<is_n; i++) {
2019:       bb[is_idx[i]] = diag*xx[is_idx[i]];
2020:     }
2021:     VecRestoreArrayRead(x,&xx);
2022:     VecRestoreArray(b,&bb);
2023:   }

2025:   /* Make a copy of the IS and  sort it */
2026:   /* allocate memory for rows,sizes */
2027:   PetscMalloc2(is_n,&rows,2*is_n,&sizes);

2029:   /* copy IS values to rows, and sort them */
2030:   for (i=0; i<is_n; i++) rows[i] = is_idx[i];
2031:   PetscSortInt(is_n,rows);

2033:   if (baij->keepnonzeropattern) {
2034:     for (i=0; i<is_n; i++) sizes[i] = 1;
2035:     bs_max          = is_n;
2036:   } else {
2037:     MatZeroRows_SeqBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);
2038:     A->nonzerostate++;
2039:   }

2041:   for (i=0,j=0; i<bs_max; j+=sizes[i],i++) {
2042:     row = rows[j];
2043:     if (row < 0 || row > A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range",row);
2044:     count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
2045:     aa    = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
2046:     if (sizes[i] == bs && !baij->keepnonzeropattern) {
2047:       if (diag != (PetscScalar)0.0) {
2048:         if (baij->ilen[row/bs] > 0) {
2049:           baij->ilen[row/bs]       = 1;
2050:           baij->j[baij->i[row/bs]] = row/bs;

2052:           PetscArrayzero(aa,count*bs);
2053:         }
2054:         /* Now insert all the diagonal values for this bs */
2055:         for (k=0; k<bs; k++) {
2056:           (*A->ops->setvalues)(A,1,rows+j+k,1,rows+j+k,&diag,INSERT_VALUES);
2057:         }
2058:       } else { /* (diag == 0.0) */
2059:         baij->ilen[row/bs] = 0;
2060:       } /* end (diag == 0.0) */
2061:     } else { /* (sizes[i] != bs) */
2062: #if defined(PETSC_USE_DEBUG)
2063:       if (sizes[i] != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal Error. Value should be 1");
2064: #endif
2065:       for (k=0; k<count; k++) {
2066:         aa[0] =  zero;
2067:         aa   += bs;
2068:       }
2069:       if (diag != (PetscScalar)0.0) {
2070:         (*A->ops->setvalues)(A,1,rows+j,1,rows+j,&diag,INSERT_VALUES);
2071:       }
2072:     }
2073:   }

2075:   PetscFree2(rows,sizes);
2076:   MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY);
2077:   return(0);
2078: }

2080: PetscErrorCode MatZeroRowsColumns_SeqBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
2081: {
2082:   Mat_SeqBAIJ       *baij=(Mat_SeqBAIJ*)A->data;
2083:   PetscErrorCode    ierr;
2084:   PetscInt          i,j,k,count;
2085:   PetscInt          bs   =A->rmap->bs,bs2=baij->bs2,row,col;
2086:   PetscScalar       zero = 0.0;
2087:   MatScalar         *aa;
2088:   const PetscScalar *xx;
2089:   PetscScalar       *bb;
2090:   PetscBool         *zeroed,vecs = PETSC_FALSE;

2093:   /* fix right hand side if needed */
2094:   if (x && b) {
2095:     VecGetArrayRead(x,&xx);
2096:     VecGetArray(b,&bb);
2097:     vecs = PETSC_TRUE;
2098:   }

2100:   /* zero the columns */
2101:   PetscCalloc1(A->rmap->n,&zeroed);
2102:   for (i=0; i<is_n; i++) {
2103:     if (is_idx[i] < 0 || is_idx[i] >= A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range",is_idx[i]);
2104:     zeroed[is_idx[i]] = PETSC_TRUE;
2105:   }
2106:   for (i=0; i<A->rmap->N; i++) {
2107:     if (!zeroed[i]) {
2108:       row = i/bs;
2109:       for (j=baij->i[row]; j<baij->i[row+1]; j++) {
2110:         for (k=0; k<bs; k++) {
2111:           col = bs*baij->j[j] + k;
2112:           if (zeroed[col]) {
2113:             aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
2114:             if (vecs) bb[i] -= aa[0]*xx[col];
2115:             aa[0] = 0.0;
2116:           }
2117:         }
2118:       }
2119:     } else if (vecs) bb[i] = diag*xx[i];
2120:   }
2121:   PetscFree(zeroed);
2122:   if (vecs) {
2123:     VecRestoreArrayRead(x,&xx);
2124:     VecRestoreArray(b,&bb);
2125:   }

2127:   /* zero the rows */
2128:   for (i=0; i<is_n; i++) {
2129:     row   = is_idx[i];
2130:     count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
2131:     aa    = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
2132:     for (k=0; k<count; k++) {
2133:       aa[0] =  zero;
2134:       aa   += bs;
2135:     }
2136:     if (diag != (PetscScalar)0.0) {
2137:       (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);
2138:     }
2139:   }
2140:   MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY);
2141:   return(0);
2142: }

2144: PetscErrorCode MatSetValues_SeqBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
2145: {
2146:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
2147:   PetscInt       *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
2148:   PetscInt       *imax=a->imax,*ai=a->i,*ailen=a->ilen;
2149:   PetscInt       *aj  =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
2151:   PetscInt       ridx,cidx,bs2=a->bs2;
2152:   PetscBool      roworiented=a->roworiented;
2153:   MatScalar      *ap=NULL,value=0.0,*aa=a->a,*bap;

2156:   for (k=0; k<m; k++) { /* loop over added rows */
2157:     row  = im[k];
2158:     brow = row/bs;
2159:     if (row < 0) continue;
2160: #if defined(PETSC_USE_DEBUG)
2161:     if (row >= A->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1);
2162: #endif
2163:     rp   = aj + ai[brow];
2164:     if (!A->structure_only) ap = aa + bs2*ai[brow];
2165:     rmax = imax[brow];
2166:     nrow = ailen[brow];
2167:     low  = 0;
2168:     high = nrow;
2169:     for (l=0; l<n; l++) { /* loop over added columns */
2170:       if (in[l] < 0) continue;
2171: #if defined(PETSC_USE_DEBUG)
2172:       if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1);
2173: #endif
2174:       col  = in[l]; bcol = col/bs;
2175:       ridx = row % bs; cidx = col % bs;
2176:       if (!A->structure_only) {
2177:         if (roworiented) {
2178:           value = v[l + k*n];
2179:         } else {
2180:           value = v[k + l*m];
2181:         }
2182:       }
2183:       if (col <= lastcol) low = 0; else high = nrow;
2184:       lastcol = col;
2185:       while (high-low > 7) {
2186:         t = (low+high)/2;
2187:         if (rp[t] > bcol) high = t;
2188:         else              low  = t;
2189:       }
2190:       for (i=low; i<high; i++) {
2191:         if (rp[i] > bcol) break;
2192:         if (rp[i] == bcol) {
2193:           bap = ap +  bs2*i + bs*cidx + ridx;
2194:           if (!A->structure_only) {
2195:             if (is == ADD_VALUES) *bap += value;
2196:             else                  *bap  = value;
2197:           }
2198:           goto noinsert1;
2199:         }
2200:       }
2201:       if (nonew == 1) goto noinsert1;
2202:       if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
2203:       if (A->structure_only) {
2204:         MatSeqXAIJReallocateAIJ_structure_only(A,a->mbs,bs2,nrow,brow,bcol,rmax,ai,aj,rp,imax,nonew,MatScalar);
2205:       } else {
2206:         MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar);
2207:       }
2208:       N = nrow++ - 1; high++;
2209:       /* shift up all the later entries in this row */
2210:       PetscArraymove(rp+i+1,rp+i,N-i+1);
2211:       rp[i] = bcol;
2212:       if (!A->structure_only) {
2213:         PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));
2214:         PetscArrayzero(ap+bs2*i,bs2);
2215:         ap[bs2*i + bs*cidx + ridx] = value;
2216:       }
2217:       a->nz++;
2218:       A->nonzerostate++;
2219: noinsert1:;
2220:       low = i;
2221:     }
2222:     ailen[brow] = nrow;
2223:   }
2224:   return(0);
2225: }

2227: PetscErrorCode MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info)
2228: {
2229:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)inA->data;
2230:   Mat            outA;
2232:   PetscBool      row_identity,col_identity;

2235:   if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels = 0 supported for in-place ILU");
2236:   ISIdentity(row,&row_identity);
2237:   ISIdentity(col,&col_identity);
2238:   if (!row_identity || !col_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Row and column permutations must be identity for in-place ILU");

2240:   outA            = inA;
2241:   inA->factortype = MAT_FACTOR_LU;
2242:   PetscFree(inA->solvertype);
2243:   PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);

2245:   MatMarkDiagonal_SeqBAIJ(inA);

2247:   PetscObjectReference((PetscObject)row);
2248:   ISDestroy(&a->row);
2249:   a->row = row;
2250:   PetscObjectReference((PetscObject)col);
2251:   ISDestroy(&a->col);
2252:   a->col = col;

2254:   /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */
2255:   ISDestroy(&a->icol);
2256:   ISInvertPermutation(col,PETSC_DECIDE,&a->icol);
2257:   PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);

2259:   MatSeqBAIJSetNumericFactorization_inplace(inA,(PetscBool)(row_identity && col_identity));
2260:   if (!a->solve_work) {
2261:     PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);
2262:     PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));
2263:   }
2264:   MatLUFactorNumeric(outA,inA,info);
2265:   return(0);
2266: }

2268: PetscErrorCode  MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,PetscInt *indices)
2269: {
2270:   Mat_SeqBAIJ *baij = (Mat_SeqBAIJ*)mat->data;
2271:   PetscInt    i,nz,mbs;

2274:   nz  = baij->maxnz;
2275:   mbs = baij->mbs;
2276:   for (i=0; i<nz; i++) {
2277:     baij->j[i] = indices[i];
2278:   }
2279:   baij->nz = nz;
2280:   for (i=0; i<mbs; i++) {
2281:     baij->ilen[i] = baij->imax[i];
2282:   }
2283:   return(0);
2284: }

2286: /*@
2287:     MatSeqBAIJSetColumnIndices - Set the column indices for all the rows
2288:        in the matrix.

2290:   Input Parameters:
2291: +  mat - the SeqBAIJ matrix
2292: -  indices - the column indices

2294:   Level: advanced

2296:   Notes:
2297:     This can be called if you have precomputed the nonzero structure of the
2298:   matrix and want to provide it to the matrix object to improve the performance
2299:   of the MatSetValues() operation.

2301:     You MUST have set the correct numbers of nonzeros per row in the call to
2302:   MatCreateSeqBAIJ(), and the columns indices MUST be sorted.

2304:     MUST be called before any calls to MatSetValues();

2306: @*/
2307: PetscErrorCode  MatSeqBAIJSetColumnIndices(Mat mat,PetscInt *indices)
2308: {

2314:   PetscUseMethod(mat,"MatSeqBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
2315:   return(0);
2316: }

2318: PetscErrorCode MatGetRowMaxAbs_SeqBAIJ(Mat A,Vec v,PetscInt idx[])
2319: {
2320:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
2322:   PetscInt       i,j,n,row,bs,*ai,*aj,mbs;
2323:   PetscReal      atmp;
2324:   PetscScalar    *x,zero = 0.0;
2325:   MatScalar      *aa;
2326:   PetscInt       ncols,brow,krow,kcol;

2329:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2330:   bs  = A->rmap->bs;
2331:   aa  = a->a;
2332:   ai  = a->i;
2333:   aj  = a->j;
2334:   mbs = a->mbs;

2336:   VecSet(v,zero);
2337:   VecGetArray(v,&x);
2338:   VecGetLocalSize(v,&n);
2339:   if (n != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector");
2340:   for (i=0; i<mbs; i++) {
2341:     ncols = ai[1] - ai[0]; ai++;
2342:     brow  = bs*i;
2343:     for (j=0; j<ncols; j++) {
2344:       for (kcol=0; kcol<bs; kcol++) {
2345:         for (krow=0; krow<bs; krow++) {
2346:           atmp = PetscAbsScalar(*aa);aa++;
2347:           row  = brow + krow;   /* row index */
2348:           if (PetscAbsScalar(x[row]) < atmp) {x[row] = atmp; if (idx) idx[row] = bs*(*aj) + kcol;}
2349:         }
2350:       }
2351:       aj++;
2352:     }
2353:   }
2354:   VecRestoreArray(v,&x);
2355:   return(0);
2356: }

2358: PetscErrorCode MatCopy_SeqBAIJ(Mat A,Mat B,MatStructure str)
2359: {

2363:   /* If the two matrices have the same copy implementation, use fast copy. */
2364:   if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
2365:     Mat_SeqBAIJ *a  = (Mat_SeqBAIJ*)A->data;
2366:     Mat_SeqBAIJ *b  = (Mat_SeqBAIJ*)B->data;
2367:     PetscInt    ambs=a->mbs,bmbs=b->mbs,abs=A->rmap->bs,bbs=B->rmap->bs,bs2=abs*abs;

2369:     if (a->i[ambs] != b->i[bmbs]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzero blocks in matrices A %D and B %D are different",a->i[ambs],b->i[bmbs]);
2370:     if (abs != bbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Block size A %D and B %D are different",abs,bbs);
2371:     PetscArraycpy(b->a,a->a,bs2*a->i[ambs]);
2372:     PetscObjectStateIncrease((PetscObject)B);
2373:   } else {
2374:     MatCopy_Basic(A,B,str);
2375:   }
2376:   return(0);
2377: }

2379: PetscErrorCode MatSetUp_SeqBAIJ(Mat A)
2380: {

2384:   MatSeqBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,0);
2385:   return(0);
2386: }

2388: static PetscErrorCode MatSeqBAIJGetArray_SeqBAIJ(Mat A,PetscScalar *array[])
2389: {
2390:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;

2393:   *array = a->a;
2394:   return(0);
2395: }

2397: static PetscErrorCode MatSeqBAIJRestoreArray_SeqBAIJ(Mat A,PetscScalar *array[])
2398: {
2400:   *array = NULL;
2401:   return(0);
2402: }

2404: PetscErrorCode MatAXPYGetPreallocation_SeqBAIJ(Mat Y,Mat X,PetscInt *nnz)
2405: {
2406:   PetscInt       bs = Y->rmap->bs,mbs = Y->rmap->N/bs;
2407:   Mat_SeqBAIJ    *x = (Mat_SeqBAIJ*)X->data;
2408:   Mat_SeqBAIJ    *y = (Mat_SeqBAIJ*)Y->data;

2412:   /* Set the number of nonzeros in the new matrix */
2413:   MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);
2414:   return(0);
2415: }

2417: PetscErrorCode MatAXPY_SeqBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
2418: {
2419:   Mat_SeqBAIJ    *x = (Mat_SeqBAIJ*)X->data,*y = (Mat_SeqBAIJ*)Y->data;
2421:   PetscInt       bs=Y->rmap->bs,bs2=bs*bs;
2422:   PetscBLASInt   one=1;

2425:   if (str == SAME_NONZERO_PATTERN) {
2426:     PetscScalar  alpha = a;
2427:     PetscBLASInt bnz;
2428:     PetscBLASIntCast(x->nz*bs2,&bnz);
2429:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one));
2430:     PetscObjectStateIncrease((PetscObject)Y);
2431:   } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
2432:     MatAXPY_Basic(Y,a,X,str);
2433:   } else {
2434:     Mat      B;
2435:     PetscInt *nnz;
2436:     if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size");
2437:     PetscMalloc1(Y->rmap->N,&nnz);
2438:     MatCreate(PetscObjectComm((PetscObject)Y),&B);
2439:     PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);
2440:     MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);
2441:     MatSetBlockSizesFromMats(B,Y,Y);
2442:     MatSetType(B,(MatType) ((PetscObject)Y)->type_name);
2443:     MatAXPYGetPreallocation_SeqBAIJ(Y,X,nnz);
2444:     MatSeqBAIJSetPreallocation(B,bs,0,nnz);
2445:     MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
2446:     MatHeaderReplace(Y,&B);
2447:     PetscFree(nnz);
2448:   }
2449:   return(0);
2450: }

2452: PetscErrorCode MatRealPart_SeqBAIJ(Mat A)
2453: {
2454:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
2455:   PetscInt    i,nz = a->bs2*a->i[a->mbs];
2456:   MatScalar   *aa = a->a;

2459:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2460:   return(0);
2461: }

2463: PetscErrorCode MatImaginaryPart_SeqBAIJ(Mat A)
2464: {
2465:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
2466:   PetscInt    i,nz = a->bs2*a->i[a->mbs];
2467:   MatScalar   *aa = a->a;

2470:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2471:   return(0);
2472: }

2474: /*
2475:     Code almost idential to MatGetColumnIJ_SeqAIJ() should share common code
2476: */
2477: PetscErrorCode MatGetColumnIJ_SeqBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool  *done)
2478: {
2479:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
2481:   PetscInt       bs = A->rmap->bs,i,*collengths,*cia,*cja,n = A->cmap->n/bs,m = A->rmap->n/bs;
2482:   PetscInt       nz = a->i[m],row,*jj,mr,col;

2485:   *nn = n;
2486:   if (!ia) return(0);
2487:   if (symmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for BAIJ matrices");
2488:   else {
2489:     PetscCalloc1(n,&collengths);
2490:     PetscMalloc1(n+1,&cia);
2491:     PetscMalloc1(nz,&cja);
2492:     jj   = a->j;
2493:     for (i=0; i<nz; i++) {
2494:       collengths[jj[i]]++;
2495:     }
2496:     cia[0] = oshift;
2497:     for (i=0; i<n; i++) {
2498:       cia[i+1] = cia[i] + collengths[i];
2499:     }
2500:     PetscArrayzero(collengths,n);
2501:     jj   = a->j;
2502:     for (row=0; row<m; row++) {
2503:       mr = a->i[row+1] - a->i[row];
2504:       for (i=0; i<mr; i++) {
2505:         col = *jj++;

2507:         cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
2508:       }
2509:     }
2510:     PetscFree(collengths);
2511:     *ia  = cia; *ja = cja;
2512:   }
2513:   return(0);
2514: }

2516: PetscErrorCode MatRestoreColumnIJ_SeqBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool  *done)
2517: {

2521:   if (!ia) return(0);
2522:   PetscFree(*ia);
2523:   PetscFree(*ja);
2524:   return(0);
2525: }

2527: /*
2528:  MatGetColumnIJ_SeqBAIJ_Color() and MatRestoreColumnIJ_SeqBAIJ_Color() are customized from
2529:  MatGetColumnIJ_SeqBAIJ() and MatRestoreColumnIJ_SeqBAIJ() by adding an output
2530:  spidx[], index of a->a, to be used in MatTransposeColoringCreate() and MatFDColoringCreate()
2531:  */
2532: PetscErrorCode MatGetColumnIJ_SeqBAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool  *done)
2533: {
2534:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
2536:   PetscInt       i,*collengths,*cia,*cja,n=a->nbs,m=a->mbs;
2537:   PetscInt       nz = a->i[m],row,*jj,mr,col;
2538:   PetscInt       *cspidx;

2541:   *nn = n;
2542:   if (!ia) return(0);

2544:   PetscCalloc1(n,&collengths);
2545:   PetscMalloc1(n+1,&cia);
2546:   PetscMalloc1(nz,&cja);
2547:   PetscMalloc1(nz,&cspidx);
2548:   jj   = a->j;
2549:   for (i=0; i<nz; i++) {
2550:     collengths[jj[i]]++;
2551:   }
2552:   cia[0] = oshift;
2553:   for (i=0; i<n; i++) {
2554:     cia[i+1] = cia[i] + collengths[i];
2555:   }
2556:   PetscArrayzero(collengths,n);
2557:   jj   = a->j;
2558:   for (row=0; row<m; row++) {
2559:     mr = a->i[row+1] - a->i[row];
2560:     for (i=0; i<mr; i++) {
2561:       col = *jj++;
2562:       cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */
2563:       cja[cia[col] + collengths[col]++ - oshift]  = row + oshift;
2564:     }
2565:   }
2566:   PetscFree(collengths);
2567:   *ia    = cia;
2568:   *ja    = cja;
2569:   *spidx = cspidx;
2570:   return(0);
2571: }

2573: PetscErrorCode MatRestoreColumnIJ_SeqBAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool  *done)
2574: {

2578:   MatRestoreColumnIJ_SeqBAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);
2579:   PetscFree(*spidx);
2580:   return(0);
2581: }

2583: PetscErrorCode MatShift_SeqBAIJ(Mat Y,PetscScalar a)
2584: {
2586:   Mat_SeqBAIJ     *aij = (Mat_SeqBAIJ*)Y->data;

2589:   if (!Y->preallocated || !aij->nz) {
2590:     MatSeqBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);
2591:   }
2592:   MatShift_Basic(Y,a);
2593:   return(0);
2594: }

2596: /* -------------------------------------------------------------------*/
2597: static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ,
2598:                                        MatGetRow_SeqBAIJ,
2599:                                        MatRestoreRow_SeqBAIJ,
2600:                                        MatMult_SeqBAIJ_N,
2601:                                /* 4*/  MatMultAdd_SeqBAIJ_N,
2602:                                        MatMultTranspose_SeqBAIJ,
2603:                                        MatMultTransposeAdd_SeqBAIJ,
2604:                                        0,
2605:                                        0,
2606:                                        0,
2607:                                /* 10*/ 0,
2608:                                        MatLUFactor_SeqBAIJ,
2609:                                        0,
2610:                                        0,
2611:                                        MatTranspose_SeqBAIJ,
2612:                                /* 15*/ MatGetInfo_SeqBAIJ,
2613:                                        MatEqual_SeqBAIJ,
2614:                                        MatGetDiagonal_SeqBAIJ,
2615:                                        MatDiagonalScale_SeqBAIJ,
2616:                                        MatNorm_SeqBAIJ,
2617:                                /* 20*/ 0,
2618:                                        MatAssemblyEnd_SeqBAIJ,
2619:                                        MatSetOption_SeqBAIJ,
2620:                                        MatZeroEntries_SeqBAIJ,
2621:                                /* 24*/ MatZeroRows_SeqBAIJ,
2622:                                        0,
2623:                                        0,
2624:                                        0,
2625:                                        0,
2626:                                /* 29*/ MatSetUp_SeqBAIJ,
2627:                                        0,
2628:                                        0,
2629:                                        0,
2630:                                        0,
2631:                                /* 34*/ MatDuplicate_SeqBAIJ,
2632:                                        0,
2633:                                        0,
2634:                                        MatILUFactor_SeqBAIJ,
2635:                                        0,
2636:                                /* 39*/ MatAXPY_SeqBAIJ,
2637:                                        MatCreateSubMatrices_SeqBAIJ,
2638:                                        MatIncreaseOverlap_SeqBAIJ,
2639:                                        MatGetValues_SeqBAIJ,
2640:                                        MatCopy_SeqBAIJ,
2641:                                /* 44*/ 0,
2642:                                        MatScale_SeqBAIJ,
2643:                                        MatShift_SeqBAIJ,
2644:                                        0,
2645:                                        MatZeroRowsColumns_SeqBAIJ,
2646:                                /* 49*/ 0,
2647:                                        MatGetRowIJ_SeqBAIJ,
2648:                                        MatRestoreRowIJ_SeqBAIJ,
2649:                                        MatGetColumnIJ_SeqBAIJ,
2650:                                        MatRestoreColumnIJ_SeqBAIJ,
2651:                                /* 54*/ MatFDColoringCreate_SeqXAIJ,
2652:                                        0,
2653:                                        0,
2654:                                        0,
2655:                                        MatSetValuesBlocked_SeqBAIJ,
2656:                                /* 59*/ MatCreateSubMatrix_SeqBAIJ,
2657:                                        MatDestroy_SeqBAIJ,
2658:                                        MatView_SeqBAIJ,
2659:                                        0,
2660:                                        0,
2661:                                /* 64*/ 0,
2662:                                        0,
2663:                                        0,
2664:                                        0,
2665:                                        0,
2666:                                /* 69*/ MatGetRowMaxAbs_SeqBAIJ,
2667:                                        0,
2668:                                        MatConvert_Basic,
2669:                                        0,
2670:                                        0,
2671:                                /* 74*/ 0,
2672:                                        MatFDColoringApply_BAIJ,
2673:                                        0,
2674:                                        0,
2675:                                        0,
2676:                                /* 79*/ 0,
2677:                                        0,
2678:                                        0,
2679:                                        0,
2680:                                        MatLoad_SeqBAIJ,
2681:                                /* 84*/ 0,
2682:                                        0,
2683:                                        0,
2684:                                        0,
2685:                                        0,
2686:                                /* 89*/ 0,
2687:                                        0,
2688:                                        0,
2689:                                        0,
2690:                                        0,
2691:                                /* 94*/ 0,
2692:                                        0,
2693:                                        0,
2694:                                        0,
2695:                                        0,
2696:                                /* 99*/ 0,
2697:                                        0,
2698:                                        0,
2699:                                        0,
2700:                                        0,
2701:                                /*104*/ 0,
2702:                                        MatRealPart_SeqBAIJ,
2703:                                        MatImaginaryPart_SeqBAIJ,
2704:                                        0,
2705:                                        0,
2706:                                /*109*/ 0,
2707:                                        0,
2708:                                        0,
2709:                                        0,
2710:                                        MatMissingDiagonal_SeqBAIJ,
2711:                                /*114*/ 0,
2712:                                        0,
2713:                                        0,
2714:                                        0,
2715:                                        0,
2716:                                /*119*/ 0,
2717:                                        0,
2718:                                        MatMultHermitianTranspose_SeqBAIJ,
2719:                                        MatMultHermitianTransposeAdd_SeqBAIJ,
2720:                                        0,
2721:                                /*124*/ 0,
2722:                                        0,
2723:                                        MatInvertBlockDiagonal_SeqBAIJ,
2724:                                        0,
2725:                                        0,
2726:                                /*129*/ 0,
2727:                                        0,
2728:                                        0,
2729:                                        0,
2730:                                        0,
2731:                                /*134*/ 0,
2732:                                        0,
2733:                                        0,
2734:                                        0,
2735:                                        0,
2736:                                /*139*/ MatSetBlockSizes_Default,
2737:                                        0,
2738:                                        0,
2739:                                        MatFDColoringSetUp_SeqXAIJ,
2740:                                        0,
2741:                                 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqBAIJ,
2742:                                        MatDestroySubMatrices_SeqBAIJ
2743: };

2745: PetscErrorCode  MatStoreValues_SeqBAIJ(Mat mat)
2746: {
2747:   Mat_SeqBAIJ    *aij = (Mat_SeqBAIJ*)mat->data;
2748:   PetscInt       nz   = aij->i[aij->mbs]*aij->bs2;

2752:   if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");

2754:   /* allocate space for values if not already there */
2755:   if (!aij->saved_values) {
2756:     PetscMalloc1(nz+1,&aij->saved_values);
2757:     PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));
2758:   }

2760:   /* copy values over */
2761:   PetscArraycpy(aij->saved_values,aij->a,nz);
2762:   return(0);
2763: }

2765: PetscErrorCode  MatRetrieveValues_SeqBAIJ(Mat mat)
2766: {
2767:   Mat_SeqBAIJ    *aij = (Mat_SeqBAIJ*)mat->data;
2769:   PetscInt       nz = aij->i[aij->mbs]*aij->bs2;

2772:   if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
2773:   if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first");

2775:   /* copy values over */
2776:   PetscArraycpy(aij->a,aij->saved_values,nz);
2777:   return(0);
2778: }

2780: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*);
2781: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat, MatType,MatReuse,Mat*);

2783: PetscErrorCode  MatSeqBAIJSetPreallocation_SeqBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
2784: {
2785:   Mat_SeqBAIJ    *b;
2787:   PetscInt       i,mbs,nbs,bs2;
2788:   PetscBool      flg = PETSC_FALSE,skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE;

2791:   if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
2792:   if (nz == MAT_SKIP_ALLOCATION) {
2793:     skipallocation = PETSC_TRUE;
2794:     nz             = 0;
2795:   }

2797:   MatSetBlockSize(B,PetscAbs(bs));
2798:   PetscLayoutSetUp(B->rmap);
2799:   PetscLayoutSetUp(B->cmap);
2800:   PetscLayoutGetBlockSize(B->rmap,&bs);

2802:   B->preallocated = PETSC_TRUE;

2804:   mbs = B->rmap->n/bs;
2805:   nbs = B->cmap->n/bs;
2806:   bs2 = bs*bs;

2808:   if (mbs*bs!=B->rmap->n || nbs*bs!=B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows %D, cols %D must be divisible by blocksize %D",B->rmap->N,B->cmap->n,bs);

2810:   if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
2811:   if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
2812:   if (nnz) {
2813:     for (i=0; i<mbs; i++) {
2814:       if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]);
2815:       if (nnz[i] > nbs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D rowlength %D",i,nnz[i],nbs);
2816:     }
2817:   }

2819:   b    = (Mat_SeqBAIJ*)B->data;
2820:   PetscOptionsBegin(PetscObjectComm((PetscObject)B),NULL,"Optimize options for SEQBAIJ matrix 2 ","Mat");
2821:   PetscOptionsBool("-mat_no_unroll","Do not optimize for block size (slow)",NULL,flg,&flg,NULL);
2822:   PetscOptionsEnd();

2824:   if (!flg) {
2825:     switch (bs) {
2826:     case 1:
2827:       B->ops->mult    = MatMult_SeqBAIJ_1;
2828:       B->ops->multadd = MatMultAdd_SeqBAIJ_1;
2829:       break;
2830:     case 2:
2831:       B->ops->mult    = MatMult_SeqBAIJ_2;
2832:       B->ops->multadd = MatMultAdd_SeqBAIJ_2;
2833:       break;
2834:     case 3:
2835:       B->ops->mult    = MatMult_SeqBAIJ_3;
2836:       B->ops->multadd = MatMultAdd_SeqBAIJ_3;
2837:       break;
2838:     case 4:
2839:       B->ops->mult    = MatMult_SeqBAIJ_4;
2840:       B->ops->multadd = MatMultAdd_SeqBAIJ_4;
2841:       break;
2842:     case 5:
2843:       B->ops->mult    = MatMult_SeqBAIJ_5;
2844:       B->ops->multadd = MatMultAdd_SeqBAIJ_5;
2845:       break;
2846:     case 6:
2847:       B->ops->mult    = MatMult_SeqBAIJ_6;
2848:       B->ops->multadd = MatMultAdd_SeqBAIJ_6;
2849:       break;
2850:     case 7:
2851:       B->ops->mult    = MatMult_SeqBAIJ_7;
2852:       B->ops->multadd = MatMultAdd_SeqBAIJ_7;
2853:       break;
2854:     case 9:
2855: #if defined(PETSC_HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && defined(PETSC_USE_REAL_DOUBLE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_64BIT_INDICES)
2856:       B->ops->mult    = MatMult_SeqBAIJ_9_AVX2;
2857:       B->ops->multadd = MatMultAdd_SeqBAIJ_9_AVX2;
2858: #else
2859:       B->ops->mult    = MatMult_SeqBAIJ_N;
2860:       B->ops->multadd = MatMultAdd_SeqBAIJ_N;
2861: #endif
2862:       break;
2863:     case 11:
2864:       B->ops->mult    = MatMult_SeqBAIJ_11;
2865:       B->ops->multadd = MatMultAdd_SeqBAIJ_11;
2866:       break;
2867:     case 15:
2868:       B->ops->mult    = MatMult_SeqBAIJ_15_ver1;
2869:       B->ops->multadd = MatMultAdd_SeqBAIJ_N;
2870:       break;
2871:     default:
2872:       B->ops->mult    = MatMult_SeqBAIJ_N;
2873:       B->ops->multadd = MatMultAdd_SeqBAIJ_N;
2874:       break;
2875:     }
2876:   }
2877:   B->ops->sor = MatSOR_SeqBAIJ;
2878:   b->mbs = mbs;
2879:   b->nbs = nbs;
2880:   if (!skipallocation) {
2881:     if (!b->imax) {
2882:       PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);
2883:       PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));

2885:       b->free_imax_ilen = PETSC_TRUE;
2886:     }
2887:     /* b->ilen will count nonzeros in each block row so far. */
2888:     for (i=0; i<mbs; i++) b->ilen[i] = 0;
2889:     if (!nnz) {
2890:       if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
2891:       else if (nz < 0) nz = 1;
2892:       nz = PetscMin(nz,nbs);
2893:       for (i=0; i<mbs; i++) b->imax[i] = nz;
2894:       nz = nz*mbs;
2895:     } else {
2896:       PetscInt64 nz64 = 0;
2897:       for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];}
2898:       PetscIntCast(nz64,&nz);
2899:     }

2901:     /* allocate the matrix space */
2902:     MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
2903:     if (B->structure_only) {
2904:       PetscMalloc1(nz,&b->j);
2905:       PetscMalloc1(B->rmap->N+1,&b->i);
2906:       PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));
2907:     } else {
2908:       PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);
2909:       PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));
2910:       PetscArrayzero(b->a,nz*bs2);
2911:     }
2912:     PetscArrayzero(b->j,nz);

2914:     if (B->structure_only) {
2915:       b->singlemalloc = PETSC_FALSE;
2916:       b->free_a       = PETSC_FALSE;
2917:     } else {
2918:       b->singlemalloc = PETSC_TRUE;
2919:       b->free_a       = PETSC_TRUE;
2920:     }
2921:     b->free_ij = PETSC_TRUE;

2923:     b->i[0] = 0;
2924:     for (i=1; i<mbs+1; i++) {
2925:       b->i[i] = b->i[i-1] + b->imax[i-1];
2926:     }

2928:   } else {
2929:     b->free_a  = PETSC_FALSE;
2930:     b->free_ij = PETSC_FALSE;
2931:   }

2933:   b->bs2              = bs2;
2934:   b->mbs              = mbs;
2935:   b->nz               = 0;
2936:   b->maxnz            = nz;
2937:   B->info.nz_unneeded = (PetscReal)b->maxnz*bs2;
2938:   B->was_assembled    = PETSC_FALSE;
2939:   B->assembled        = PETSC_FALSE;
2940:   if (realalloc) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);}
2941:   return(0);
2942: }

2944: PetscErrorCode MatSeqBAIJSetPreallocationCSR_SeqBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[],const PetscScalar V[])
2945: {
2946:   PetscInt       i,m,nz,nz_max=0,*nnz;
2947:   PetscScalar    *values=0;
2948:   PetscBool      roworiented = ((Mat_SeqBAIJ*)B->data)->roworiented;

2952:   if (bs < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs);
2953:   PetscLayoutSetBlockSize(B->rmap,bs);
2954:   PetscLayoutSetBlockSize(B->cmap,bs);
2955:   PetscLayoutSetUp(B->rmap);
2956:   PetscLayoutSetUp(B->cmap);
2957:   PetscLayoutGetBlockSize(B->rmap,&bs);
2958:   m    = B->rmap->n/bs;

2960:   if (ii[0] != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %D",ii[0]);
2961:   PetscMalloc1(m+1, &nnz);
2962:   for (i=0; i<m; i++) {
2963:     nz = ii[i+1]- ii[i];
2964:     if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D",i,nz);
2965:     nz_max = PetscMax(nz_max, nz);
2966:     nnz[i] = nz;
2967:   }
2968:   MatSeqBAIJSetPreallocation(B,bs,0,nnz);
2969:   PetscFree(nnz);

2971:   values = (PetscScalar*)V;
2972:   if (!values) {
2973:     PetscCalloc1(bs*bs*(nz_max+1),&values);
2974:   }
2975:   for (i=0; i<m; i++) {
2976:     PetscInt          ncols  = ii[i+1] - ii[i];
2977:     const PetscInt    *icols = jj + ii[i];
2978:     if (bs == 1 || !roworiented) {
2979:       const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0);
2980:       MatSetValuesBlocked_SeqBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);
2981:     } else {
2982:       PetscInt j;
2983:       for (j=0; j<ncols; j++) {
2984:         const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0);
2985:         MatSetValuesBlocked_SeqBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);
2986:       }
2987:     }
2988:   }
2989:   if (!V) { PetscFree(values); }
2990:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2991:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2992:   MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
2993:   return(0);
2994: }

2996: /*@C
2997:    MatSeqBAIJGetArray - gives access to the array where the data for a MATSEQBAIJ matrix is stored

2999:    Not Collective

3001:    Input Parameter:
3002: .  mat - a MATSEQBAIJ matrix

3004:    Output Parameter:
3005: .   array - pointer to the data

3007:    Level: intermediate

3009: .seealso: MatSeqBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
3010: @*/
3011: PetscErrorCode MatSeqBAIJGetArray(Mat A,PetscScalar **array)
3012: {

3016:   PetscUseMethod(A,"MatSeqBAIJGetArray_C",(Mat,PetscScalar**),(A,array));
3017:   return(0);
3018: }

3020: /*@C
3021:    MatSeqBAIJRestoreArray - returns access to the array where the data for a MATSEQBAIJ matrix is stored obtained by MatSeqBAIJGetArray()

3023:    Not Collective

3025:    Input Parameters:
3026: +  mat - a MATSEQBAIJ matrix
3027: -  array - pointer to the data

3029:    Level: intermediate

3031: .seealso: MatSeqBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray()
3032: @*/
3033: PetscErrorCode MatSeqBAIJRestoreArray(Mat A,PetscScalar **array)
3034: {

3038:   PetscUseMethod(A,"MatSeqBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));
3039:   return(0);
3040: }

3042: /*MC
3043:    MATSEQBAIJ - MATSEQBAIJ = "seqbaij" - A matrix type to be used for sequential block sparse matrices, based on
3044:    block sparse compressed row format.

3046:    Options Database Keys:
3047: . -mat_type seqbaij - sets the matrix type to "seqbaij" during a call to MatSetFromOptions()

3049:    Level: beginner

3051:    Notes:
3052:     MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no
3053:     space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored

3055: .seealso: MatCreateSeqBAIJ()
3056: M*/

3058: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqBSTRM(Mat, MatType,MatReuse,Mat*);

3060: PETSC_EXTERN PetscErrorCode MatCreate_SeqBAIJ(Mat B)
3061: {
3063:   PetscMPIInt    size;
3064:   Mat_SeqBAIJ    *b;

3067:   MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);
3068:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1");

3070:   PetscNewLog(B,&b);
3071:   B->data = (void*)b;
3072:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));

3074:   b->row          = 0;
3075:   b->col          = 0;
3076:   b->icol         = 0;
3077:   b->reallocs     = 0;
3078:   b->saved_values = 0;

3080:   b->roworiented        = PETSC_TRUE;
3081:   b->nonew              = 0;
3082:   b->diag               = 0;
3083:   B->spptr              = 0;
3084:   B->info.nz_unneeded   = (PetscReal)b->maxnz*b->bs2;
3085:   b->keepnonzeropattern = PETSC_FALSE;

3087:   PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJGetArray_C",MatSeqBAIJGetArray_SeqBAIJ);
3088:   PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJRestoreArray_C",MatSeqBAIJRestoreArray_SeqBAIJ);
3089:   PetscObjectComposeFunction((PetscObject)B,"MatInvertBlockDiagonal_C",MatInvertBlockDiagonal_SeqBAIJ);
3090:   PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqBAIJ);
3091:   PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqBAIJ);
3092:   PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetColumnIndices_C",MatSeqBAIJSetColumnIndices_SeqBAIJ);
3093:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaij_seqaij_C",MatConvert_SeqBAIJ_SeqAIJ);
3094:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaij_seqsbaij_C",MatConvert_SeqBAIJ_SeqSBAIJ);
3095:   PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetPreallocation_C",MatSeqBAIJSetPreallocation_SeqBAIJ);
3096:   PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetPreallocationCSR_C",MatSeqBAIJSetPreallocationCSR_SeqBAIJ);
3097:   PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqBAIJ);
3098: #if defined(PETSC_HAVE_HYPRE)
3099:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaij_hypre_C",MatConvert_AIJ_HYPRE);
3100: #endif
3101:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaij_is_C",MatConvert_XAIJ_IS);
3102:   PetscObjectChangeTypeName((PetscObject)B,MATSEQBAIJ);
3103:   return(0);
3104: }

3106: PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace)
3107: {
3108:   Mat_SeqBAIJ    *c = (Mat_SeqBAIJ*)C->data,*a = (Mat_SeqBAIJ*)A->data;
3110:   PetscInt       i,mbs = a->mbs,nz = a->nz,bs2 = a->bs2;

3113:   if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix");

3115:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3116:     c->imax           = a->imax;
3117:     c->ilen           = a->ilen;
3118:     c->free_imax_ilen = PETSC_FALSE;
3119:   } else {
3120:     PetscMalloc2(mbs,&c->imax,mbs,&c->ilen);
3121:     PetscLogObjectMemory((PetscObject)C,2*mbs*sizeof(PetscInt));
3122:     for (i=0; i<mbs; i++) {
3123:       c->imax[i] = a->imax[i];
3124:       c->ilen[i] = a->ilen[i];
3125:     }
3126:     c->free_imax_ilen = PETSC_TRUE;
3127:   }

3129:   /* allocate the matrix space */
3130:   if (mallocmatspace) {
3131:     if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3132:       PetscCalloc1(bs2*nz,&c->a);
3133:       PetscLogObjectMemory((PetscObject)C,a->i[mbs]*bs2*sizeof(PetscScalar));

3135:       c->i            = a->i;
3136:       c->j            = a->j;
3137:       c->singlemalloc = PETSC_FALSE;
3138:       c->free_a       = PETSC_TRUE;
3139:       c->free_ij      = PETSC_FALSE;
3140:       c->parent       = A;
3141:       C->preallocated = PETSC_TRUE;
3142:       C->assembled    = PETSC_TRUE;

3144:       PetscObjectReference((PetscObject)A);
3145:       MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
3146:       MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
3147:     } else {
3148:       PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);
3149:       PetscLogObjectMemory((PetscObject)C,a->i[mbs]*(bs2*sizeof(PetscScalar)+sizeof(PetscInt))+(mbs+1)*sizeof(PetscInt));

3151:       c->singlemalloc = PETSC_TRUE;
3152:       c->free_a       = PETSC_TRUE;
3153:       c->free_ij      = PETSC_TRUE;

3155:       PetscArraycpy(c->i,a->i,mbs+1);
3156:       if (mbs > 0) {
3157:         PetscArraycpy(c->j,a->j,nz);
3158:         if (cpvalues == MAT_COPY_VALUES) {
3159:           PetscArraycpy(c->a,a->a,bs2*nz);
3160:         } else {
3161:           PetscArrayzero(c->a,bs2*nz);
3162:         }
3163:       }
3164:       C->preallocated = PETSC_TRUE;
3165:       C->assembled    = PETSC_TRUE;
3166:     }
3167:   }

3169:   c->roworiented = a->roworiented;
3170:   c->nonew       = a->nonew;

3172:   PetscLayoutReference(A->rmap,&C->rmap);
3173:   PetscLayoutReference(A->cmap,&C->cmap);

3175:   c->bs2         = a->bs2;
3176:   c->mbs         = a->mbs;
3177:   c->nbs         = a->nbs;

3179:   if (a->diag) {
3180:     if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3181:       c->diag      = a->diag;
3182:       c->free_diag = PETSC_FALSE;
3183:     } else {
3184:       PetscMalloc1(mbs+1,&c->diag);
3185:       PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt));
3186:       for (i=0; i<mbs; i++) c->diag[i] = a->diag[i];
3187:       c->free_diag = PETSC_TRUE;
3188:     }
3189:   } else c->diag = 0;

3191:   c->nz         = a->nz;
3192:   c->maxnz      = a->nz;         /* Since we allocate exactly the right amount */
3193:   c->solve_work = NULL;
3194:   c->mult_work  = NULL;
3195:   c->sor_workt  = NULL;
3196:   c->sor_work   = NULL;

3198:   c->compressedrow.use   = a->compressedrow.use;
3199:   c->compressedrow.nrows = a->compressedrow.nrows;
3200:   if (a->compressedrow.use) {
3201:     i    = a->compressedrow.nrows;
3202:     PetscMalloc2(i+1,&c->compressedrow.i,i+1,&c->compressedrow.rindex);
3203:     PetscLogObjectMemory((PetscObject)C,(2*i+1)*sizeof(PetscInt));
3204:     PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1);
3205:     PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i);
3206:   } else {
3207:     c->compressedrow.use    = PETSC_FALSE;
3208:     c->compressedrow.i      = NULL;
3209:     c->compressedrow.rindex = NULL;
3210:   }
3211:   C->nonzerostate = A->nonzerostate;

3213:   PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
3214:   return(0);
3215: }

3217: PetscErrorCode MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
3218: {

3222:   MatCreate(PetscObjectComm((PetscObject)A),B);
3223:   MatSetSizes(*B,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);
3224:   MatSetType(*B,MATSEQBAIJ);
3225:   MatDuplicateNoCreate_SeqBAIJ(*B,A,cpvalues,PETSC_TRUE);
3226:   return(0);
3227: }

3229: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
3230: PetscErrorCode MatLoad_SeqBAIJ_Binary(Mat mat,PetscViewer viewer)
3231: {
3232:   PetscInt       header[4],M,N,nz,bs,m,n,mbs,nbs,rows,cols,sum,i,j,k;
3233:   PetscInt       *rowidxs,*colidxs;
3234:   PetscScalar    *matvals;

3238:   PetscViewerSetUp(viewer);

3240:   /* read matrix header */
3241:   PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);
3242:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file");
3243:   M = header[1]; N = header[2]; nz = header[3];
3244:   if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M);
3245:   if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N);
3246:   if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqBAIJ");

3248:   /* set block sizes from the viewer's .info file */
3249:   MatLoad_Binary_BlockSizes(mat,viewer);
3250:   /* set local and global sizes if not set already */
3251:   if (mat->rmap->n < 0) mat->rmap->n = M;
3252:   if (mat->cmap->n < 0) mat->cmap->n = N;
3253:   if (mat->rmap->N < 0) mat->rmap->N = M;
3254:   if (mat->cmap->N < 0) mat->cmap->N = N;
3255:   PetscLayoutSetUp(mat->rmap);
3256:   PetscLayoutSetUp(mat->cmap);

3258:   /* check if the matrix sizes are correct */
3259:   MatGetSize(mat,&rows,&cols);
3260:   if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%D, %D) than the input matrix (%D, %D)",M,N,rows,cols);
3261:   MatGetBlockSize(mat,&bs);
3262:   MatGetLocalSize(mat,&m,&n);
3263:   mbs = m/bs; nbs = n/bs;

3265:   /* read in row lengths, column indices and nonzero values */
3266:   PetscMalloc1(m+1,&rowidxs);
3267:   PetscViewerBinaryRead(viewer,rowidxs+1,m,NULL,PETSC_INT);
3268:   rowidxs[0] = 0; for (i=0; i<m; i++) rowidxs[i+1] += rowidxs[i];
3269:   sum = rowidxs[m];
3270:   if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent matrix data in file: nonzeros = %D, sum-row-lengths = %D\n",nz,sum);

3272:   /* read in column indices and nonzero values */
3273:   PetscMalloc2(rowidxs[m],&colidxs,nz,&matvals);
3274:   PetscViewerBinaryRead(viewer,colidxs,rowidxs[m],NULL,PETSC_INT);
3275:   PetscViewerBinaryRead(viewer,matvals,rowidxs[m],NULL,PETSC_SCALAR);

3277:   { /* preallocate matrix storage */
3278:     PetscBT   bt; /* helper bit set to count nonzeros */
3279:     PetscInt  *nnz;
3280:     PetscBool sbaij;

3282:     PetscBTCreate(nbs,&bt);
3283:     PetscCalloc1(mbs,&nnz);
3284:     PetscObjectTypeCompare((PetscObject)mat,MATSEQSBAIJ,&sbaij);
3285:     for (i=0; i<mbs; i++) {
3286:       PetscBTMemzero(nbs,bt);
3287:       for (k=0; k<bs; k++) {
3288:         PetscInt row = bs*i + k;
3289:         for (j=rowidxs[row]; j<rowidxs[row+1]; j++) {
3290:           PetscInt col = colidxs[j];
3291:           if (!sbaij || col >= row)
3292:             if (!PetscBTLookupSet(bt,col/bs)) nnz[i]++;
3293:         }
3294:       }
3295:     }
3296:     PetscBTDestroy(&bt);
3297:     MatSeqBAIJSetPreallocation(mat,bs,0,nnz);
3298:     MatSeqSBAIJSetPreallocation(mat,bs,0,nnz);
3299:     PetscFree(nnz);
3300:   }

3302:   /* store matrix values */
3303:   for (i=0; i<m; i++) {
3304:     PetscInt row = i, s = rowidxs[i], e = rowidxs[i+1];
3305:     (*mat->ops->setvalues)(mat,1,&row,e-s,colidxs+s,matvals+s,INSERT_VALUES);
3306:   }

3308:   PetscFree(rowidxs);
3309:   PetscFree2(colidxs,matvals);
3310:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
3311:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
3312:   return(0);
3313: }

3315: PetscErrorCode MatLoad_SeqBAIJ(Mat mat,PetscViewer viewer)
3316: {
3318:   PetscBool      isbinary;

3321:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3322:   if (!isbinary) SETERRQ2(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)mat)->type_name);
3323:   MatLoad_SeqBAIJ_Binary(mat,viewer);
3324:   return(0);
3325: }

3327: /*@C
3328:    MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block
3329:    compressed row) format.  For good matrix assembly performance the
3330:    user should preallocate the matrix storage by setting the parameter nz
3331:    (or the array nnz).  By setting these parameters accurately, performance
3332:    during matrix assembly can be increased by more than a factor of 50.

3334:    Collective

3336:    Input Parameters:
3337: +  comm - MPI communicator, set to PETSC_COMM_SELF
3338: .  bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
3339:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
3340: .  m - number of rows
3341: .  n - number of columns
3342: .  nz - number of nonzero blocks  per block row (same for all rows)
3343: -  nnz - array containing the number of nonzero blocks in the various block rows
3344:          (possibly different for each block row) or NULL

3346:    Output Parameter:
3347: .  A - the matrix

3349:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
3350:    MatXXXXSetPreallocation() paradigm instead of this routine directly.
3351:    [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]

3353:    Options Database Keys:
3354: +   -mat_no_unroll - uses code that does not unroll the loops in the
3355:                      block calculations (much slower)
3356: -    -mat_block_size - size of the blocks to use

3358:    Level: intermediate

3360:    Notes:
3361:    The number of rows and columns must be divisible by blocksize.

3363:    If the nnz parameter is given then the nz parameter is ignored

3365:    A nonzero block is any block that as 1 or more nonzeros in it

3367:    The block AIJ format is fully compatible with standard Fortran 77
3368:    storage.  That is, the stored row and column indices can begin at
3369:    either one (as in Fortran) or zero.  See the users' manual for details.

3371:    Specify the preallocated storage with either nz or nnz (not both).
3372:    Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3373:    allocation.  See Users-Manual: ch_mat for details.
3374:    matrices.

3376: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateBAIJ()
3377: @*/
3378: PetscErrorCode  MatCreateSeqBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
3379: {

3383:   MatCreate(comm,A);
3384:   MatSetSizes(*A,m,n,m,n);
3385:   MatSetType(*A,MATSEQBAIJ);
3386:   MatSeqBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);
3387:   return(0);
3388: }

3390: /*@C
3391:    MatSeqBAIJSetPreallocation - Sets the block size and expected nonzeros
3392:    per row in the matrix. For good matrix assembly performance the
3393:    user should preallocate the matrix storage by setting the parameter nz
3394:    (or the array nnz).  By setting these parameters accurately, performance
3395:    during matrix assembly can be increased by more than a factor of 50.

3397:    Collective

3399:    Input Parameters:
3400: +  B - the matrix
3401: .  bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
3402:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
3403: .  nz - number of block nonzeros per block row (same for all rows)
3404: -  nnz - array containing the number of block nonzeros in the various block rows
3405:          (possibly different for each block row) or NULL

3407:    Options Database Keys:
3408: +   -mat_no_unroll - uses code that does not unroll the loops in the
3409:                      block calculations (much slower)
3410: -   -mat_block_size - size of the blocks to use

3412:    Level: intermediate

3414:    Notes:
3415:    If the nnz parameter is given then the nz parameter is ignored

3417:    You can call MatGetInfo() to get information on how effective the preallocation was;
3418:    for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3419:    You can also run with the option -info and look for messages with the string
3420:    malloc in them to see if additional memory allocation was needed.

3422:    The block AIJ format is fully compatible with standard Fortran 77
3423:    storage.  That is, the stored row and column indices can begin at
3424:    either one (as in Fortran) or zero.  See the users' manual for details.

3426:    Specify the preallocated storage with either nz or nnz (not both).
3427:    Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
3428:    allocation.  See Users-Manual: ch_mat for details.

3430: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateBAIJ(), MatGetInfo()
3431: @*/
3432: PetscErrorCode  MatSeqBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
3433: {

3440:   PetscTryMethod(B,"MatSeqBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
3441:   return(0);
3442: }

3444: /*@C
3445:    MatSeqBAIJSetPreallocationCSR - Creates a sparse parallel matrix in BAIJ format using the given nonzero structure and (optional) numerical values

3447:    Collective

3449:    Input Parameters:
3450: +  B - the matrix
3451: .  i - the indices into j for the start of each local row (starts with zero)
3452: .  j - the column indices for each local row (starts with zero) these must be sorted for each row
3453: -  v - optional values in the matrix

3455:    Level: advanced

3457:    Notes:
3458:    The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED.  For example, C programs
3459:    may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is
3460:    over rows within a block and the last index is over columns within a block row.  Fortran programs will likely set
3461:    MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
3462:    block column and the second index is over columns within a block.

3464:    Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries and usually the numerical values as well

3466: .seealso: MatCreate(), MatCreateSeqBAIJ(), MatSetValues(), MatSeqBAIJSetPreallocation(), MATSEQBAIJ
3467: @*/
3468: PetscErrorCode  MatSeqBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[])
3469: {

3476:   PetscTryMethod(B,"MatSeqBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));
3477:   return(0);
3478: }


3481: /*@
3482:      MatCreateSeqBAIJWithArrays - Creates an sequential BAIJ matrix using matrix elements provided by the user.

3484:      Collective

3486:    Input Parameters:
3487: +  comm - must be an MPI communicator of size 1
3488: .  bs - size of block
3489: .  m - number of rows
3490: .  n - number of columns
3491: .  i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row block row of the matrix
3492: .  j - column indices
3493: -  a - matrix values

3495:    Output Parameter:
3496: .  mat - the matrix

3498:    Level: advanced

3500:    Notes:
3501:        The i, j, and a arrays are not copied by this routine, the user must free these arrays
3502:     once the matrix is destroyed

3504:        You cannot set new nonzero locations into this matrix, that will generate an error.

3506:        The i and j indices are 0 based

3508:        When block size is greater than 1 the matrix values must be stored using the BAIJ storage format (see the BAIJ code to determine this).

3510:       The order of the entries in values is the same as the block compressed sparse row storage format; that is, it is
3511:       the same as a three dimensional array in Fortran values(bs,bs,nnz) that contains the first column of the first
3512:       block, followed by the second column of the first block etc etc.  That is, the blocks are contiguous in memory
3513:       with column-major ordering within blocks.

3515: .seealso: MatCreate(), MatCreateBAIJ(), MatCreateSeqBAIJ()

3517: @*/
3518: PetscErrorCode  MatCreateSeqBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
3519: {
3521:   PetscInt       ii;
3522:   Mat_SeqBAIJ    *baij;

3525:   if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs);
3526:   if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0");

3528:   MatCreate(comm,mat);
3529:   MatSetSizes(*mat,m,n,m,n);
3530:   MatSetType(*mat,MATSEQBAIJ);
3531:   MatSeqBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,0);
3532:   baij = (Mat_SeqBAIJ*)(*mat)->data;
3533:   PetscMalloc2(m,&baij->imax,m,&baij->ilen);
3534:   PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));

3536:   baij->i = i;
3537:   baij->j = j;
3538:   baij->a = a;

3540:   baij->singlemalloc = PETSC_FALSE;
3541:   baij->nonew        = -1;             /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
3542:   baij->free_a       = PETSC_FALSE;
3543:   baij->free_ij      = PETSC_FALSE;

3545:   for (ii=0; ii<m; ii++) {
3546:     baij->ilen[ii] = baij->imax[ii] = i[ii+1] - i[ii];
3547: #if defined(PETSC_USE_DEBUG)
3548:     if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]);
3549: #endif
3550:   }
3551: #if defined(PETSC_USE_DEBUG)
3552:   for (ii=0; ii<baij->i[m]; ii++) {
3553:     if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
3554:     if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]);
3555:   }
3556: #endif

3558:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
3559:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
3560:   return(0);
3561: }

3563: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
3564: {
3566:   PetscMPIInt    size;

3569:   MPI_Comm_size(comm,&size);
3570:   if (size == 1 && scall == MAT_REUSE_MATRIX) {
3571:     MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
3572:   } else {
3573:     MatCreateMPIMatConcatenateSeqMat_MPIBAIJ(comm,inmat,n,scall,outmat);
3574:   }
3575:   return(0);
3576: }