Actual source code: baij.c

petsc-master 2019-09-15
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*);
 19: PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*);

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

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

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

 52:       if (PetscAbsScalar(diag[0] + shift) < PETSC_MACHINE_EPSILON) {
 53:         if (allowzeropivot) {
 54:           A->factorerrortype             = MAT_FACTOR_NUMERIC_ZEROPIVOT;
 55:           A->factorerror_zeropivot_value = PetscAbsScalar(diag[0]);
 56:           A->factorerror_zeropivot_row   = i;
 57:           PetscInfo1(A,"Zero pivot, row %D\n",i);
 58:         } 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);
 59:       }

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

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

148:   its = its*lits;
149:   if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat");
150:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
151:   if (fshift) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Sorry, no support for diagonal shift");
152:   if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Sorry, no support for non-trivial relaxation factor");
153:   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");

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

1213:   PetscObjectChangeTypeName((PetscObject)A,0);
1214:   PetscObjectComposeFunction((PetscObject)A,"MatInvertBlockDiagonal_C",NULL);
1215:   PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);
1216:   PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);
1217:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBAIJSetColumnIndices_C",NULL);
1218:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_seqaij_C",NULL);
1219:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_seqsbaij_C",NULL);
1220:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",NULL);
1221:   PetscObjectComposeFunction((PetscObject)A,"MatSeqBAIJSetPreallocationCSR_C",NULL);
1222:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_seqbstrm_C",NULL);
1223:   PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);
1224: #if defined(PETSC_HAVE_HYPRE)
1225:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_hypre_C",NULL);
1226: #endif
1227:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqbaij_is_C",NULL);
1228:   PetscObjectComposeFunction((PetscObject)A,"MatPtAP_is_seqaij_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: static PetscErrorCode MatView_SeqBAIJ_Binary(Mat A,PetscViewer viewer)
1418: {
1419:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
1421:   PetscInt       i,*col_lens,bs = A->rmap->bs,count,*jj,j,k,l,bs2=a->bs2;
1422:   int            fd;
1423:   PetscScalar    *aa;
1424:   FILE           *file;

1427:   PetscViewerBinaryGetDescriptor(viewer,&fd);
1428:   PetscMalloc1(4+A->rmap->N,&col_lens);
1429:   col_lens[0] = MAT_FILE_CLASSID;

1431:   col_lens[1] = A->rmap->N;
1432:   col_lens[2] = A->cmap->n;
1433:   col_lens[3] = a->nz*bs2;

1435:   /* store lengths of each row and write (including header) to file */
1436:   count = 0;
1437:   for (i=0; i<a->mbs; i++) {
1438:     for (j=0; j<bs; j++) {
1439:       col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]);
1440:     }
1441:   }
1442:   PetscBinaryWrite(fd,col_lens,4+A->rmap->N,PETSC_INT,PETSC_TRUE);
1443:   PetscFree(col_lens);

1445:   /* store column indices (zero start index) */
1446:   PetscMalloc1((a->nz+1)*bs2,&jj);
1447:   count = 0;
1448:   for (i=0; i<a->mbs; i++) {
1449:     for (j=0; j<bs; j++) {
1450:       for (k=a->i[i]; k<a->i[i+1]; k++) {
1451:         for (l=0; l<bs; l++) {
1452:           jj[count++] = bs*a->j[k] + l;
1453:         }
1454:       }
1455:     }
1456:   }
1457:   PetscBinaryWrite(fd,jj,bs2*a->nz,PETSC_INT,PETSC_FALSE);
1458:   PetscFree(jj);

1460:   /* store nonzero values */
1461:   PetscMalloc1((a->nz+1)*bs2,&aa);
1462:   count = 0;
1463:   for (i=0; i<a->mbs; i++) {
1464:     for (j=0; j<bs; j++) {
1465:       for (k=a->i[i]; k<a->i[i+1]; k++) {
1466:         for (l=0; l<bs; l++) {
1467:           aa[count++] = a->a[bs2*k + l*bs + j];
1468:         }
1469:       }
1470:     }
1471:   }
1472:   PetscBinaryWrite(fd,aa,bs2*a->nz,PETSC_SCALAR,PETSC_FALSE);
1473:   PetscFree(aa);

1475:   PetscViewerBinaryGetInfoPointer(viewer,&file);
1476:   if (file) {
1477:     fprintf(file,"-matload_block_size %d\n",(int)A->rmap->bs);
1478:   }
1479:   return(0);
1480: }

1482: static PetscErrorCode MatView_SeqBAIJ_ASCII_structonly(Mat A,PetscViewer viewer)
1483: {
1485:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
1486:   PetscInt       i,bs = A->rmap->bs,k;

1489:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1490:   for (i=0; i<a->mbs; i++) {
1491:     PetscViewerASCIIPrintf(viewer,"row %D-%D:",i*bs,i*bs+bs-1);
1492:     for (k=a->i[i]; k<a->i[i+1]; k++) {
1493:       PetscViewerASCIIPrintf(viewer," (%D-%D) ",bs*a->j[k],bs*a->j[k]+bs-1);
1494:     }
1495:     PetscViewerASCIIPrintf(viewer,"\n");
1496:   }
1497:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1498:   return(0);
1499: }

1501: static PetscErrorCode MatView_SeqBAIJ_ASCII(Mat A,PetscViewer viewer)
1502: {
1503:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
1504:   PetscErrorCode    ierr;
1505:   PetscInt          i,j,bs = A->rmap->bs,k,l,bs2=a->bs2;
1506:   PetscViewerFormat format;

1509:   if (A->structure_only) {
1510:     MatView_SeqBAIJ_ASCII_structonly(A,viewer);
1511:     return(0);
1512:   }

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

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

1599:   PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);
1600:   PetscViewerGetFormat(viewer,&format);
1601:   PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);

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

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

1657:     for (i=0; i<a->nz*a->bs2; i++) {
1658:       if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]);
1659:     }
1660:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1661:     PetscDrawGetPopup(draw,&popup);
1662:     PetscDrawScalePopup(popup,0.0,maxv);

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

1684: static PetscErrorCode MatView_SeqBAIJ_Draw(Mat A,PetscViewer viewer)
1685: {
1687:   PetscReal      xl,yl,xr,yr,w,h;
1688:   PetscDraw      draw;
1689:   PetscBool      isnull;

1692:   PetscViewerDrawGetDraw(viewer,0,&draw);
1693:   PetscDrawIsNull(draw,&isnull);
1694:   if (isnull) return(0);

1696:   xr   = A->cmap->n; yr = A->rmap->N; h = yr/10.0; w = xr/10.0;
1697:   xr  += w;          yr += h;        xl = -w;     yl = -h;
1698:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
1699:   PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);
1700:   PetscDrawZoom(draw,MatView_SeqBAIJ_Draw_Zoom,A);
1701:   PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);
1702:   PetscDrawSave(draw);
1703:   return(0);
1704: }

1706: PetscErrorCode MatView_SeqBAIJ(Mat A,PetscViewer viewer)
1707: {
1709:   PetscBool      iascii,isbinary,isdraw;

1712:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1713:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1714:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1715:   if (iascii) {
1716:     MatView_SeqBAIJ_ASCII(A,viewer);
1717:   } else if (isbinary) {
1718:     MatView_SeqBAIJ_Binary(A,viewer);
1719:   } else if (isdraw) {
1720:     MatView_SeqBAIJ_Draw(A,viewer);
1721:   } else {
1722:     Mat B;
1723:     MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
1724:     MatView(B,viewer);
1725:     MatDestroy(&B);
1726:   }
1727:   return(0);
1728: }


1731: PetscErrorCode MatGetValues_SeqBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[])
1732: {
1733:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
1734:   PetscInt    *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j;
1735:   PetscInt    *ai = a->i,*ailen = a->ilen;
1736:   PetscInt    brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2;
1737:   MatScalar   *ap,*aa = a->a;

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

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

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

1899: PetscErrorCode MatAssemblyEnd_SeqBAIJ(Mat A,MatAssemblyType mode)
1900: {
1901:   Mat_SeqBAIJ    *a     = (Mat_SeqBAIJ*)A->data;
1902:   PetscInt       fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax;
1903:   PetscInt       m      = A->rmap->N,*ip,N,*ailen = a->ilen;
1905:   PetscInt       mbs  = a->mbs,bs2 = a->bs2,rmax = 0;
1906:   MatScalar      *aa  = a->a,*ap;
1907:   PetscReal      ratio=0.6;

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

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

1933:   /* reset ilen and imax for each row */
1934:   a->nonzerorowcnt = 0;
1935:   if (A->structure_only) {
1936:     PetscFree2(a->imax,a->ilen);
1937:   } else { /* !A->structure_only */
1938:     for (i=0; i<mbs; i++) {
1939:       ailen[i] = imax[i] = ai[i+1] - ai[i];
1940:       a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0);
1941:     }
1942:   }
1943:   a->nz = ai[mbs];

1945:   /* diagonals may have moved, so kill the diagonal pointers */
1946:   a->idiagvalid = PETSC_FALSE;
1947:   if (fshift && a->diag) {
1948:     PetscFree(a->diag);
1949:     PetscLogObjectMemory((PetscObject)A,-(mbs+1)*sizeof(PetscInt));
1950:     a->diag = 0;
1951:   }
1952:   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);
1953:   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);
1954:   PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);
1955:   PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);

1957:   A->info.mallocs    += a->reallocs;
1958:   a->reallocs         = 0;
1959:   A->info.nz_unneeded = (PetscReal)fshift*bs2;
1960:   a->rmax             = rmax;

1962:   if (!A->structure_only) {
1963:     MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,mbs,ratio);
1964:   }
1965:   return(0);
1966: }

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

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

2009: PetscErrorCode MatZeroRows_SeqBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
2010: {
2011:   Mat_SeqBAIJ       *baij=(Mat_SeqBAIJ*)A->data;
2012:   PetscErrorCode    ierr;
2013:   PetscInt          i,j,k,count,*rows;
2014:   PetscInt          bs=A->rmap->bs,bs2=baij->bs2,*sizes,row,bs_max;
2015:   PetscScalar       zero = 0.0;
2016:   MatScalar         *aa;
2017:   const PetscScalar *xx;
2018:   PetscScalar       *bb;

2021:   /* fix right hand side if needed */
2022:   if (x && b) {
2023:     VecGetArrayRead(x,&xx);
2024:     VecGetArray(b,&bb);
2025:     for (i=0; i<is_n; i++) {
2026:       bb[is_idx[i]] = diag*xx[is_idx[i]];
2027:     }
2028:     VecRestoreArrayRead(x,&xx);
2029:     VecRestoreArray(b,&bb);
2030:   }

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

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

2040:   if (baij->keepnonzeropattern) {
2041:     for (i=0; i<is_n; i++) sizes[i] = 1;
2042:     bs_max          = is_n;
2043:   } else {
2044:     MatZeroRows_SeqBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);
2045:     A->nonzerostate++;
2046:   }

2048:   for (i=0,j=0; i<bs_max; j+=sizes[i],i++) {
2049:     row = rows[j];
2050:     if (row < 0 || row > A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range",row);
2051:     count = (baij->i[row/bs +1] - baij->i[row/bs])*bs;
2052:     aa    = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs);
2053:     if (sizes[i] == bs && !baij->keepnonzeropattern) {
2054:       if (diag != (PetscScalar)0.0) {
2055:         if (baij->ilen[row/bs] > 0) {
2056:           baij->ilen[row/bs]       = 1;
2057:           baij->j[baij->i[row/bs]] = row/bs;

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

2082:   PetscFree2(rows,sizes);
2083:   MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY);
2084:   return(0);
2085: }

2087: PetscErrorCode MatZeroRowsColumns_SeqBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b)
2088: {
2089:   Mat_SeqBAIJ       *baij=(Mat_SeqBAIJ*)A->data;
2090:   PetscErrorCode    ierr;
2091:   PetscInt          i,j,k,count;
2092:   PetscInt          bs   =A->rmap->bs,bs2=baij->bs2,row,col;
2093:   PetscScalar       zero = 0.0;
2094:   MatScalar         *aa;
2095:   const PetscScalar *xx;
2096:   PetscScalar       *bb;
2097:   PetscBool         *zeroed,vecs = PETSC_FALSE;

2100:   /* fix right hand side if needed */
2101:   if (x && b) {
2102:     VecGetArrayRead(x,&xx);
2103:     VecGetArray(b,&bb);
2104:     vecs = PETSC_TRUE;
2105:   }

2107:   /* zero the columns */
2108:   PetscCalloc1(A->rmap->n,&zeroed);
2109:   for (i=0; i<is_n; i++) {
2110:     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]);
2111:     zeroed[is_idx[i]] = PETSC_TRUE;
2112:   }
2113:   for (i=0; i<A->rmap->N; i++) {
2114:     if (!zeroed[i]) {
2115:       row = i/bs;
2116:       for (j=baij->i[row]; j<baij->i[row+1]; j++) {
2117:         for (k=0; k<bs; k++) {
2118:           col = bs*baij->j[j] + k;
2119:           if (zeroed[col]) {
2120:             aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k;
2121:             if (vecs) bb[i] -= aa[0]*xx[col];
2122:             aa[0] = 0.0;
2123:           }
2124:         }
2125:       }
2126:     } else if (vecs) bb[i] = diag*xx[i];
2127:   }
2128:   PetscFree(zeroed);
2129:   if (vecs) {
2130:     VecRestoreArrayRead(x,&xx);
2131:     VecRestoreArray(b,&bb);
2132:   }

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

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

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

2234: PetscErrorCode MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info)
2235: {
2236:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)inA->data;
2237:   Mat            outA;
2239:   PetscBool      row_identity,col_identity;

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

2247:   outA            = inA;
2248:   inA->factortype = MAT_FACTOR_LU;
2249:   PetscFree(inA->solvertype);
2250:   PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);

2252:   MatMarkDiagonal_SeqBAIJ(inA);

2254:   PetscObjectReference((PetscObject)row);
2255:   ISDestroy(&a->row);
2256:   a->row = row;
2257:   PetscObjectReference((PetscObject)col);
2258:   ISDestroy(&a->col);
2259:   a->col = col;

2261:   /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */
2262:   ISDestroy(&a->icol);
2263:   ISInvertPermutation(col,PETSC_DECIDE,&a->icol);
2264:   PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);

2266:   MatSeqBAIJSetNumericFactorization_inplace(inA,(PetscBool)(row_identity && col_identity));
2267:   if (!a->solve_work) {
2268:     PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);
2269:     PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));
2270:   }
2271:   MatLUFactorNumeric(outA,inA,info);
2272:   return(0);
2273: }

2275: PetscErrorCode  MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,PetscInt *indices)
2276: {
2277:   Mat_SeqBAIJ *baij = (Mat_SeqBAIJ*)mat->data;
2278:   PetscInt    i,nz,mbs;

2281:   nz  = baij->maxnz;
2282:   mbs = baij->mbs;
2283:   for (i=0; i<nz; i++) {
2284:     baij->j[i] = indices[i];
2285:   }
2286:   baij->nz = nz;
2287:   for (i=0; i<mbs; i++) {
2288:     baij->ilen[i] = baij->imax[i];
2289:   }
2290:   return(0);
2291: }

2293: /*@
2294:     MatSeqBAIJSetColumnIndices - Set the column indices for all the rows
2295:        in the matrix.

2297:   Input Parameters:
2298: +  mat - the SeqBAIJ matrix
2299: -  indices - the column indices

2301:   Level: advanced

2303:   Notes:
2304:     This can be called if you have precomputed the nonzero structure of the
2305:   matrix and want to provide it to the matrix object to improve the performance
2306:   of the MatSetValues() operation.

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

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

2313: @*/
2314: PetscErrorCode  MatSeqBAIJSetColumnIndices(Mat mat,PetscInt *indices)
2315: {

2321:   PetscUseMethod(mat,"MatSeqBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));
2322:   return(0);
2323: }

2325: PetscErrorCode MatGetRowMaxAbs_SeqBAIJ(Mat A,Vec v,PetscInt idx[])
2326: {
2327:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
2329:   PetscInt       i,j,n,row,bs,*ai,*aj,mbs;
2330:   PetscReal      atmp;
2331:   PetscScalar    *x,zero = 0.0;
2332:   MatScalar      *aa;
2333:   PetscInt       ncols,brow,krow,kcol;

2336:   if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2337:   bs  = A->rmap->bs;
2338:   aa  = a->a;
2339:   ai  = a->i;
2340:   aj  = a->j;
2341:   mbs = a->mbs;

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

2365: PetscErrorCode MatCopy_SeqBAIJ(Mat A,Mat B,MatStructure str)
2366: {

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

2376:     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]);
2377:     if (abs != bbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Block size A %D and B %D are different",abs,bbs);
2378:     PetscArraycpy(b->a,a->a,bs2*a->i[ambs]);
2379:     PetscObjectStateIncrease((PetscObject)B);
2380:   } else {
2381:     MatCopy_Basic(A,B,str);
2382:   }
2383:   return(0);
2384: }

2386: PetscErrorCode MatSetUp_SeqBAIJ(Mat A)
2387: {

2391:   MatSeqBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,0);
2392:   return(0);
2393: }

2395: PetscErrorCode MatSeqBAIJGetArray_SeqBAIJ(Mat A,PetscScalar *array[])
2396: {
2397:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;

2400:   *array = a->a;
2401:   return(0);
2402: }

2404: PetscErrorCode MatSeqBAIJRestoreArray_SeqBAIJ(Mat A,PetscScalar *array[])
2405: {
2407:   return(0);
2408: }

2410: PetscErrorCode MatAXPYGetPreallocation_SeqBAIJ(Mat Y,Mat X,PetscInt *nnz)
2411: {
2412:   PetscInt       bs = Y->rmap->bs,mbs = Y->rmap->N/bs;
2413:   Mat_SeqBAIJ    *x = (Mat_SeqBAIJ*)X->data;
2414:   Mat_SeqBAIJ    *y = (Mat_SeqBAIJ*)Y->data;

2418:   /* Set the number of nonzeros in the new matrix */
2419:   MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);
2420:   return(0);
2421: }

2423: PetscErrorCode MatAXPY_SeqBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str)
2424: {
2425:   Mat_SeqBAIJ    *x = (Mat_SeqBAIJ*)X->data,*y = (Mat_SeqBAIJ*)Y->data;
2427:   PetscInt       bs=Y->rmap->bs,bs2=bs*bs;
2428:   PetscBLASInt   one=1;

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

2458: PetscErrorCode MatRealPart_SeqBAIJ(Mat A)
2459: {
2460:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
2461:   PetscInt    i,nz = a->bs2*a->i[a->mbs];
2462:   MatScalar   *aa = a->a;

2465:   for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]);
2466:   return(0);
2467: }

2469: PetscErrorCode MatImaginaryPart_SeqBAIJ(Mat A)
2470: {
2471:   Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
2472:   PetscInt    i,nz = a->bs2*a->i[a->mbs];
2473:   MatScalar   *aa = a->a;

2476:   for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2477:   return(0);
2478: }

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

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

2513:         cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
2514:       }
2515:     }
2516:     PetscFree(collengths);
2517:     *ia  = cia; *ja = cja;
2518:   }
2519:   return(0);
2520: }

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

2527:   if (!ia) return(0);
2528:   PetscFree(*ia);
2529:   PetscFree(*ja);
2530:   return(0);
2531: }

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

2547:   *nn = n;
2548:   if (!ia) return(0);

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

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

2584:   MatRestoreColumnIJ_SeqBAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done);
2585:   PetscFree(*spidx);
2586:   return(0);
2587: }

2589: PetscErrorCode MatShift_SeqBAIJ(Mat Y,PetscScalar a)
2590: {
2592:   Mat_SeqBAIJ     *aij = (Mat_SeqBAIJ*)Y->data;

2595:   if (!Y->preallocated || !aij->nz) {
2596:     MatSeqBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);
2597:   }
2598:   MatShift_Basic(Y,a);
2599:   return(0);
2600: }

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

2751: PetscErrorCode  MatStoreValues_SeqBAIJ(Mat mat)
2752: {
2753:   Mat_SeqBAIJ    *aij = (Mat_SeqBAIJ*)mat->data;
2754:   PetscInt       nz   = aij->i[aij->mbs]*aij->bs2;

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

2760:   /* allocate space for values if not already there */
2761:   if (!aij->saved_values) {
2762:     PetscMalloc1(nz+1,&aij->saved_values);
2763:     PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));
2764:   }

2766:   /* copy values over */
2767:   PetscArraycpy(aij->saved_values,aij->a,nz);
2768:   return(0);
2769: }

2771: PetscErrorCode  MatRetrieveValues_SeqBAIJ(Mat mat)
2772: {
2773:   Mat_SeqBAIJ    *aij = (Mat_SeqBAIJ*)mat->data;
2775:   PetscInt       nz = aij->i[aij->mbs]*aij->bs2;

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

2781:   /* copy values over */
2782:   PetscArraycpy(aij->a,aij->saved_values,nz);
2783:   return(0);
2784: }

2786: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*);
2787: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat, MatType,MatReuse,Mat*);

2789: PetscErrorCode  MatSeqBAIJSetPreallocation_SeqBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz)
2790: {
2791:   Mat_SeqBAIJ    *b;
2793:   PetscInt       i,mbs,nbs,bs2;
2794:   PetscBool      flg = PETSC_FALSE,skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE;

2797:   if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
2798:   if (nz == MAT_SKIP_ALLOCATION) {
2799:     skipallocation = PETSC_TRUE;
2800:     nz             = 0;
2801:   }

2803:   MatSetBlockSize(B,PetscAbs(bs));
2804:   PetscLayoutSetUp(B->rmap);
2805:   PetscLayoutSetUp(B->cmap);
2806:   PetscLayoutGetBlockSize(B->rmap,&bs);

2808:   B->preallocated = PETSC_TRUE;

2810:   mbs = B->rmap->n/bs;
2811:   nbs = B->cmap->n/bs;
2812:   bs2 = bs*bs;

2814:   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);

2816:   if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
2817:   if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz);
2818:   if (nnz) {
2819:     for (i=0; i<mbs; i++) {
2820:       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]);
2821:       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);
2822:     }
2823:   }

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

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

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

2907:     /* allocate the matrix space */
2908:     MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);
2909:     if (B->structure_only) {
2910:       PetscMalloc1(nz,&b->j);
2911:       PetscMalloc1(B->rmap->N+1,&b->i);
2912:       PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));
2913:     } else {
2914:       PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);
2915:       PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));
2916:       PetscArrayzero(b->a,nz*bs2);
2917:     }
2918:     PetscArrayzero(b->j,nz);

2920:     if (B->structure_only) {
2921:       b->singlemalloc = PETSC_FALSE;
2922:       b->free_a       = PETSC_FALSE;
2923:     } else {
2924:       b->singlemalloc = PETSC_TRUE;
2925:       b->free_a       = PETSC_TRUE;
2926:     }
2927:     b->free_ij = PETSC_TRUE;

2929:     b->i[0] = 0;
2930:     for (i=1; i<mbs+1; i++) {
2931:       b->i[i] = b->i[i-1] + b->imax[i-1];
2932:     }

2934:   } else {
2935:     b->free_a  = PETSC_FALSE;
2936:     b->free_ij = PETSC_FALSE;
2937:   }

2939:   b->bs2              = bs2;
2940:   b->mbs              = mbs;
2941:   b->nz               = 0;
2942:   b->maxnz            = nz;
2943:   B->info.nz_unneeded = (PetscReal)b->maxnz*bs2;
2944:   B->was_assembled    = PETSC_FALSE;
2945:   B->assembled        = PETSC_FALSE;
2946:   if (realalloc) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);}
2947:   return(0);
2948: }

2950: PetscErrorCode MatSeqBAIJSetPreallocationCSR_SeqBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[],const PetscScalar V[])
2951: {
2952:   PetscInt       i,m,nz,nz_max=0,*nnz;
2953:   PetscScalar    *values=0;
2954:   PetscBool      roworiented = ((Mat_SeqBAIJ*)B->data)->roworiented;

2958:   if (bs < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs);
2959:   PetscLayoutSetBlockSize(B->rmap,bs);
2960:   PetscLayoutSetBlockSize(B->cmap,bs);
2961:   PetscLayoutSetUp(B->rmap);
2962:   PetscLayoutSetUp(B->cmap);
2963:   PetscLayoutGetBlockSize(B->rmap,&bs);
2964:   m    = B->rmap->n/bs;

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

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

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

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

3009:    Level: beginner

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

3015: .seealso: MatCreateSeqBAIJ()
3016: M*/

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

3020: PETSC_EXTERN PetscErrorCode MatCreate_SeqBAIJ(Mat B)
3021: {
3023:   PetscMPIInt    size;
3024:   Mat_SeqBAIJ    *b;

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

3030:   PetscNewLog(B,&b);
3031:   B->data = (void*)b;
3032:   PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));

3034:   b->row          = 0;
3035:   b->col          = 0;
3036:   b->icol         = 0;
3037:   b->reallocs     = 0;
3038:   b->saved_values = 0;

3040:   b->roworiented        = PETSC_TRUE;
3041:   b->nonew              = 0;
3042:   b->diag               = 0;
3043:   B->spptr              = 0;
3044:   B->info.nz_unneeded   = (PetscReal)b->maxnz*b->bs2;
3045:   b->keepnonzeropattern = PETSC_FALSE;

3047:   PetscObjectComposeFunction((PetscObject)B,"MatInvertBlockDiagonal_C",MatInvertBlockDiagonal_SeqBAIJ);
3048:   PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqBAIJ);
3049:   PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqBAIJ);
3050:   PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetColumnIndices_C",MatSeqBAIJSetColumnIndices_SeqBAIJ);
3051:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaij_seqaij_C",MatConvert_SeqBAIJ_SeqAIJ);
3052:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaij_seqsbaij_C",MatConvert_SeqBAIJ_SeqSBAIJ);
3053:   PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetPreallocation_C",MatSeqBAIJSetPreallocation_SeqBAIJ);
3054:   PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetPreallocationCSR_C",MatSeqBAIJSetPreallocationCSR_SeqBAIJ);
3055:   PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqBAIJ);
3056: #if defined(PETSC_HAVE_HYPRE)
3057:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaij_hypre_C",MatConvert_AIJ_HYPRE);
3058: #endif
3059:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqbaij_is_C",MatConvert_XAIJ_IS);
3060:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_seqbaij_C",MatPtAP_IS_XAIJ);
3061:   PetscObjectChangeTypeName((PetscObject)B,MATSEQBAIJ);
3062:   return(0);
3063: }

3065: PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace)
3066: {
3067:   Mat_SeqBAIJ    *c = (Mat_SeqBAIJ*)C->data,*a = (Mat_SeqBAIJ*)A->data;
3069:   PetscInt       i,mbs = a->mbs,nz = a->nz,bs2 = a->bs2;

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

3074:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3075:     c->imax           = a->imax;
3076:     c->ilen           = a->ilen;
3077:     c->free_imax_ilen = PETSC_FALSE;
3078:   } else {
3079:     PetscMalloc2(mbs,&c->imax,mbs,&c->ilen);
3080:     PetscLogObjectMemory((PetscObject)C,2*mbs*sizeof(PetscInt));
3081:     for (i=0; i<mbs; i++) {
3082:       c->imax[i] = a->imax[i];
3083:       c->ilen[i] = a->ilen[i];
3084:     }
3085:     c->free_imax_ilen = PETSC_TRUE;
3086:   }

3088:   /* allocate the matrix space */
3089:   if (mallocmatspace) {
3090:     if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3091:       PetscCalloc1(bs2*nz,&c->a);
3092:       PetscLogObjectMemory((PetscObject)C,a->i[mbs]*bs2*sizeof(PetscScalar));

3094:       c->i            = a->i;
3095:       c->j            = a->j;
3096:       c->singlemalloc = PETSC_FALSE;
3097:       c->free_a       = PETSC_TRUE;
3098:       c->free_ij      = PETSC_FALSE;
3099:       c->parent       = A;
3100:       C->preallocated = PETSC_TRUE;
3101:       C->assembled    = PETSC_TRUE;

3103:       PetscObjectReference((PetscObject)A);
3104:       MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
3105:       MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
3106:     } else {
3107:       PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);
3108:       PetscLogObjectMemory((PetscObject)C,a->i[mbs]*(bs2*sizeof(PetscScalar)+sizeof(PetscInt))+(mbs+1)*sizeof(PetscInt));

3110:       c->singlemalloc = PETSC_TRUE;
3111:       c->free_a       = PETSC_TRUE;
3112:       c->free_ij      = PETSC_TRUE;

3114:       PetscArraycpy(c->i,a->i,mbs+1);
3115:       if (mbs > 0) {
3116:         PetscArraycpy(c->j,a->j,nz);
3117:         if (cpvalues == MAT_COPY_VALUES) {
3118:           PetscArraycpy(c->a,a->a,bs2*nz);
3119:         } else {
3120:           PetscArrayzero(c->a,bs2*nz);
3121:         }
3122:       }
3123:       C->preallocated = PETSC_TRUE;
3124:       C->assembled    = PETSC_TRUE;
3125:     }
3126:   }

3128:   c->roworiented = a->roworiented;
3129:   c->nonew       = a->nonew;

3131:   PetscLayoutReference(A->rmap,&C->rmap);
3132:   PetscLayoutReference(A->cmap,&C->cmap);

3134:   c->bs2         = a->bs2;
3135:   c->mbs         = a->mbs;
3136:   c->nbs         = a->nbs;

3138:   if (a->diag) {
3139:     if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3140:       c->diag      = a->diag;
3141:       c->free_diag = PETSC_FALSE;
3142:     } else {
3143:       PetscMalloc1(mbs+1,&c->diag);
3144:       PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt));
3145:       for (i=0; i<mbs; i++) c->diag[i] = a->diag[i];
3146:       c->free_diag = PETSC_TRUE;
3147:     }
3148:   } else c->diag = 0;

3150:   c->nz         = a->nz;
3151:   c->maxnz      = a->nz;         /* Since we allocate exactly the right amount */
3152:   c->solve_work = NULL;
3153:   c->mult_work  = NULL;
3154:   c->sor_workt  = NULL;
3155:   c->sor_work   = NULL;

3157:   c->compressedrow.use   = a->compressedrow.use;
3158:   c->compressedrow.nrows = a->compressedrow.nrows;
3159:   if (a->compressedrow.use) {
3160:     i    = a->compressedrow.nrows;
3161:     PetscMalloc2(i+1,&c->compressedrow.i,i+1,&c->compressedrow.rindex);
3162:     PetscLogObjectMemory((PetscObject)C,(2*i+1)*sizeof(PetscInt));
3163:     PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1);
3164:     PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i);
3165:   } else {
3166:     c->compressedrow.use    = PETSC_FALSE;
3167:     c->compressedrow.i      = NULL;
3168:     c->compressedrow.rindex = NULL;
3169:   }
3170:   C->nonzerostate = A->nonzerostate;

3172:   PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);
3173:   return(0);
3174: }

3176: PetscErrorCode MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B)
3177: {

3181:   MatCreate(PetscObjectComm((PetscObject)A),B);
3182:   MatSetSizes(*B,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);
3183:   MatSetType(*B,MATSEQBAIJ);
3184:   MatDuplicateNoCreate_SeqBAIJ(*B,A,cpvalues,PETSC_TRUE);
3185:   return(0);
3186: }

3188: PetscErrorCode MatLoad_SeqBAIJ(Mat newmat,PetscViewer viewer)
3189: {
3190:   Mat_SeqBAIJ    *a;
3192:   PetscInt       i,nz,header[4],*rowlengths=0,M,N,bs = newmat->rmap->bs;
3193:   PetscInt       *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount;
3194:   PetscInt       kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols;
3195:   PetscInt       *masked,nmask,tmp,bs2,ishift;
3196:   PetscMPIInt    size;
3197:   int            fd;
3198:   PetscScalar    *aa;
3199:   MPI_Comm       comm;
3200:   PetscBool      isbinary;

3203:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3204:   if (!isbinary) SETERRQ2(PetscObjectComm((PetscObject)newmat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newmat)->type_name);

3206:   /* force binary viewer to load .info file if it has not yet done so */
3207:   PetscViewerSetUp(viewer);
3208:   PetscObjectGetComm((PetscObject)viewer,&comm);
3209:   PetscOptionsBegin(comm,NULL,"Options for loading SEQBAIJ matrix","Mat");
3210:   PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);
3211:   PetscOptionsEnd();
3212:   if (bs < 0) bs = 1;
3213:   bs2  = bs*bs;

3215:   MPI_Comm_size(comm,&size);
3216:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor");
3217:   PetscViewerBinaryGetDescriptor(viewer,&fd);
3218:   PetscBinaryRead(fd,header,4,NULL,PETSC_INT);
3219:   if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object");
3220:   M = header[1]; N = header[2]; nz = header[3];

3222:   if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqBAIJ");
3223:   if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices");

3225:   /*
3226:      This code adds extra rows to make sure the number of rows is
3227:     divisible by the blocksize
3228:   */
3229:   mbs        = M/bs;
3230:   extra_rows = bs - M + bs*(mbs);
3231:   if (extra_rows == bs) extra_rows = 0;
3232:   else mbs++;
3233:   if (extra_rows) {
3234:     PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");
3235:   }

3237:   /* Set global sizes if not already set */
3238:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
3239:     MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);
3240:   } else { /* Check if the matrix global sizes are correct */
3241:     MatGetSize(newmat,&rows,&cols);
3242:     if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */
3243:       MatGetLocalSize(newmat,&rows,&cols);
3244:     }
3245:     if (M != rows ||  N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols);
3246:   }

3248:   /* read in row lengths */
3249:   PetscMalloc1(M+extra_rows,&rowlengths);
3250:   PetscBinaryRead(fd,rowlengths,M,NULL,PETSC_INT);
3251:   for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1;

3253:   /* read in column indices */
3254:   PetscMalloc1(nz+extra_rows,&jj);
3255:   PetscBinaryRead(fd,jj,nz,NULL,PETSC_INT);
3256:   for (i=0; i<extra_rows; i++) jj[nz+i] = M+i;

3258:   /* loop over row lengths determining block row lengths */
3259:   PetscCalloc1(mbs,&browlengths);
3260:   PetscMalloc2(mbs,&mask,mbs,&masked);
3261:   PetscArrayzero(mask,mbs);
3262:   rowcount = 0;
3263:   nzcount  = 0;
3264:   for (i=0; i<mbs; i++) {
3265:     nmask = 0;
3266:     for (j=0; j<bs; j++) {
3267:       kmax = rowlengths[rowcount];
3268:       for (k=0; k<kmax; k++) {
3269:         tmp = jj[nzcount++]/bs;
3270:         if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;}
3271:       }
3272:       rowcount++;
3273:     }
3274:     browlengths[i] += nmask;
3275:     /* zero out the mask elements we set */
3276:     for (j=0; j<nmask; j++) mask[masked[j]] = 0;
3277:   }

3279:   /* Do preallocation  */
3280:   MatSeqBAIJSetPreallocation(newmat,bs,0,browlengths);
3281:   a    = (Mat_SeqBAIJ*)newmat->data;

3283:   /* set matrix "i" values */
3284:   a->i[0] = 0;
3285:   for (i=1; i<= mbs; i++) {
3286:     a->i[i]      = a->i[i-1] + browlengths[i-1];
3287:     a->ilen[i-1] = browlengths[i-1];
3288:   }
3289:   a->nz = 0;
3290:   for (i=0; i<mbs; i++) a->nz += browlengths[i];

3292:   /* read in nonzero values */
3293:   PetscMalloc1(nz+extra_rows,&aa);
3294:   PetscBinaryRead(fd,aa,nz,NULL,PETSC_SCALAR);
3295:   for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0;

3297:   /* set "a" and "j" values into matrix */
3298:   nzcount = 0; jcount = 0;
3299:   for (i=0; i<mbs; i++) {
3300:     nzcountb = nzcount;
3301:     nmask    = 0;
3302:     for (j=0; j<bs; j++) {
3303:       kmax = rowlengths[i*bs+j];
3304:       for (k=0; k<kmax; k++) {
3305:         tmp = jj[nzcount++]/bs;
3306:         if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;}
3307:       }
3308:     }
3309:     /* sort the masked values */
3310:     PetscSortInt(nmask,masked);

3312:     /* set "j" values into matrix */
3313:     maskcount = 1;
3314:     for (j=0; j<nmask; j++) {
3315:       a->j[jcount++]  = masked[j];
3316:       mask[masked[j]] = maskcount++;
3317:     }
3318:     /* set "a" values into matrix */
3319:     ishift = bs2*a->i[i];
3320:     for (j=0; j<bs; j++) {
3321:       kmax = rowlengths[i*bs+j];
3322:       for (k=0; k<kmax; k++) {
3323:         tmp       = jj[nzcountb]/bs;
3324:         block     = mask[tmp] - 1;
3325:         point     = jj[nzcountb] - bs*tmp;
3326:         idx       = ishift + bs2*block + j + bs*point;
3327:         a->a[idx] = (MatScalar)aa[nzcountb++];
3328:       }
3329:     }
3330:     /* zero out the mask elements we set */
3331:     for (j=0; j<nmask; j++) mask[masked[j]] = 0;
3332:   }
3333:   if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix");

3335:   PetscFree(rowlengths);
3336:   PetscFree(browlengths);
3337:   PetscFree(aa);
3338:   PetscFree(jj);
3339:   PetscFree2(mask,masked);

3341:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
3342:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
3343:   return(0);
3344: }

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

3353:    Collective

3355:    Input Parameters:
3356: +  comm - MPI communicator, set to PETSC_COMM_SELF
3357: .  bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row
3358:           blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs()
3359: .  m - number of rows
3360: .  n - number of columns
3361: .  nz - number of nonzero blocks  per block row (same for all rows)
3362: -  nnz - array containing the number of nonzero blocks in the various block rows
3363:          (possibly different for each block row) or NULL

3365:    Output Parameter:
3366: .  A - the matrix

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

3372:    Options Database Keys:
3373: +   -mat_no_unroll - uses code that does not unroll the loops in the
3374:                      block calculations (much slower)
3375: -    -mat_block_size - size of the blocks to use

3377:    Level: intermediate

3379:    Notes:
3380:    The number of rows and columns must be divisible by blocksize.

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

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

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

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

3395: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateBAIJ()
3396: @*/
3397: PetscErrorCode  MatCreateSeqBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
3398: {

3402:   MatCreate(comm,A);
3403:   MatSetSizes(*A,m,n,m,n);
3404:   MatSetType(*A,MATSEQBAIJ);
3405:   MatSeqBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);
3406:   return(0);
3407: }

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

3416:    Collective

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

3426:    Options Database Keys:
3427: +   -mat_no_unroll - uses code that does not unroll the loops in the
3428:                      block calculations (much slower)
3429: -   -mat_block_size - size of the blocks to use

3431:    Level: intermediate

3433:    Notes:
3434:    If the nnz parameter is given then the nz parameter is ignored

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

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

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

3449: .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateBAIJ(), MatGetInfo()
3450: @*/
3451: PetscErrorCode  MatSeqBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
3452: {

3459:   PetscTryMethod(B,"MatSeqBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
3460:   return(0);
3461: }

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

3466:    Collective

3468:    Input Parameters:
3469: +  B - the matrix
3470: .  i - the indices into j for the start of each local row (starts with zero)
3471: .  j - the column indices for each local row (starts with zero) these must be sorted for each row
3472: -  v - optional values in the matrix

3474:    Level: advanced

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

3483:    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

3485: .seealso: MatCreate(), MatCreateSeqBAIJ(), MatSetValues(), MatSeqBAIJSetPreallocation(), MATSEQBAIJ
3486: @*/
3487: PetscErrorCode  MatSeqBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[])
3488: {

3495:   PetscTryMethod(B,"MatSeqBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));
3496:   return(0);
3497: }


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

3503:      Collective

3505:    Input Parameters:
3506: +  comm - must be an MPI communicator of size 1
3507: .  bs - size of block
3508: .  m - number of rows
3509: .  n - number of columns
3510: .  i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row block row of the matrix
3511: .  j - column indices
3512: -  a - matrix values

3514:    Output Parameter:
3515: .  mat - the matrix

3517:    Level: advanced

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

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

3525:        The i and j indices are 0 based

3527:        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).

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

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

3536: @*/
3537: PetscErrorCode  MatCreateSeqBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat)
3538: {
3540:   PetscInt       ii;
3541:   Mat_SeqBAIJ    *baij;

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

3547:   MatCreate(comm,mat);
3548:   MatSetSizes(*mat,m,n,m,n);
3549:   MatSetType(*mat,MATSEQBAIJ);
3550:   MatSeqBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,0);
3551:   baij = (Mat_SeqBAIJ*)(*mat)->data;
3552:   PetscMalloc2(m,&baij->imax,m,&baij->ilen);
3553:   PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));

3555:   baij->i = i;
3556:   baij->j = j;
3557:   baij->a = a;

3559:   baij->singlemalloc = PETSC_FALSE;
3560:   baij->nonew        = -1;             /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
3561:   baij->free_a       = PETSC_FALSE;
3562:   baij->free_ij      = PETSC_FALSE;

3564:   for (ii=0; ii<m; ii++) {
3565:     baij->ilen[ii] = baij->imax[ii] = i[ii+1] - i[ii];
3566: #if defined(PETSC_USE_DEBUG)
3567:     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]);
3568: #endif
3569:   }
3570: #if defined(PETSC_USE_DEBUG)
3571:   for (ii=0; ii<baij->i[m]; ii++) {
3572:     if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]);
3573:     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]);
3574:   }
3575: #endif

3577:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
3578:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
3579:   return(0);
3580: }

3582: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat)
3583: {
3585:   PetscMPIInt    size;

3588:   MPI_Comm_size(comm,&size);
3589:   if (size == 1 && scall == MAT_REUSE_MATRIX) {
3590:     MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);
3591:   } else {
3592:     MatCreateMPIMatConcatenateSeqMat_MPIBAIJ(comm,inmat,n,scall,outmat);
3593:   }
3594:   return(0);
3595: }