Actual source code: inode.c
petsc-3.10.5 2019-03-28
2: /*
3: This file provides high performance routines for the Inode format (compressed sparse row)
4: by taking advantage of rows with identical nonzero structure (I-nodes).
5: */
6: #include <../src/mat/impls/aij/seq/aij.h>
8: static PetscErrorCode MatCreateColInode_Private(Mat A,PetscInt *size,PetscInt **ns)
9: {
10: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
12: PetscInt i,count,m,n,min_mn,*ns_row,*ns_col;
15: n = A->cmap->n;
16: m = A->rmap->n;
17: ns_row = a->inode.size;
19: min_mn = (m < n) ? m : n;
20: if (!ns) {
21: for (count=0,i=0; count<min_mn; count+=ns_row[i],i++) ;
22: for (; count+1 < n; count++,i++) ;
23: if (count < n) {
24: i++;
25: }
26: *size = i;
27: return(0);
28: }
29: PetscMalloc1(n+1,&ns_col);
31: /* Use the same row structure wherever feasible. */
32: for (count=0,i=0; count<min_mn; count+=ns_row[i],i++) {
33: ns_col[i] = ns_row[i];
34: }
36: /* if m < n; pad up the remainder with inode_limit */
37: for (; count+1 < n; count++,i++) {
38: ns_col[i] = 1;
39: }
40: /* The last node is the odd ball. padd it up with the remaining rows; */
41: if (count < n) {
42: ns_col[i] = n - count;
43: i++;
44: } else if (count > n) {
45: /* Adjust for the over estimation */
46: ns_col[i-1] += n - count;
47: }
48: *size = i;
49: *ns = ns_col;
50: return(0);
51: }
54: /*
55: This builds symmetric version of nonzero structure,
56: */
57: static PetscErrorCode MatGetRowIJ_SeqAIJ_Inode_Symmetric(Mat A,const PetscInt *iia[],const PetscInt *jja[],PetscInt ishift,PetscInt oshift)
58: {
59: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
61: PetscInt *work,*ia,*ja,nz,nslim_row,nslim_col,m,row,col,n;
62: PetscInt *tns,*tvc,*ns_row = a->inode.size,*ns_col,nsz,i1,i2;
63: const PetscInt *j,*jmax,*ai= a->i,*aj = a->j;
66: nslim_row = a->inode.node_count;
67: m = A->rmap->n;
68: n = A->cmap->n;
69: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatGetRowIJ_SeqAIJ_Inode_Symmetric: Matrix should be square");
71: /* Use the row_inode as column_inode */
72: nslim_col = nslim_row;
73: ns_col = ns_row;
75: /* allocate space for reformated inode structure */
76: PetscMalloc1(nslim_col+1,&tns);
77: PetscMalloc1(n+1,&tvc);
78: for (i1=0,tns[0]=0; i1<nslim_col; ++i1) tns[i1+1] = tns[i1]+ ns_row[i1];
80: for (i1=0,col=0; i1<nslim_col; ++i1) {
81: nsz = ns_col[i1];
82: for (i2=0; i2<nsz; ++i2,++col) tvc[col] = i1;
83: }
84: /* allocate space for row pointers */
85: PetscMalloc1(nslim_row+1,&ia);
86: *iia = ia;
87: PetscMemzero(ia,(nslim_row+1)*sizeof(PetscInt));
88: PetscMalloc1(nslim_row+1,&work);
90: /* determine the number of columns in each row */
91: ia[0] = oshift;
92: for (i1=0,row=0; i1<nslim_row; row+=ns_row[i1],i1++) {
94: j = aj + ai[row] + ishift;
95: jmax = aj + ai[row+1] + ishift;
96: if (j==jmax) continue; /* empty row */
97: col = *j++ + ishift;
98: i2 = tvc[col];
99: while (i2<i1 && j<jmax) { /* 1.[-xx-d-xx--] 2.[-xx-------],off-diagonal elemets */
100: ia[i1+1]++;
101: ia[i2+1]++;
102: i2++; /* Start col of next node */
103: while ((j<jmax) && ((col=*j+ishift)<tns[i2])) ++j;
104: i2 = tvc[col];
105: }
106: if (i2 == i1) ia[i2+1]++; /* now the diagonal element */
107: }
109: /* shift ia[i] to point to next row */
110: for (i1=1; i1<nslim_row+1; i1++) {
111: row = ia[i1-1];
112: ia[i1] += row;
113: work[i1-1] = row - oshift;
114: }
116: /* allocate space for column pointers */
117: nz = ia[nslim_row] + (!ishift);
118: PetscMalloc1(nz,&ja);
119: *jja = ja;
121: /* loop over lower triangular part putting into ja */
122: for (i1=0,row=0; i1<nslim_row; row += ns_row[i1],i1++) {
123: j = aj + ai[row] + ishift;
124: jmax = aj + ai[row+1] + ishift;
125: if (j==jmax) continue; /* empty row */
126: col = *j++ + ishift;
127: i2 = tvc[col];
128: while (i2<i1 && j<jmax) {
129: ja[work[i2]++] = i1 + oshift;
130: ja[work[i1]++] = i2 + oshift;
131: ++i2;
132: while ((j<jmax) && ((col=*j+ishift)< tns[i2])) ++j; /* Skip rest col indices in this node */
133: i2 = tvc[col];
134: }
135: if (i2 == i1) ja[work[i1]++] = i2 + oshift;
137: }
138: PetscFree(work);
139: PetscFree(tns);
140: PetscFree(tvc);
141: return(0);
142: }
144: /*
145: This builds nonsymmetric version of nonzero structure,
146: */
147: static PetscErrorCode MatGetRowIJ_SeqAIJ_Inode_Nonsymmetric(Mat A,const PetscInt *iia[],const PetscInt *jja[],PetscInt ishift,PetscInt oshift)
148: {
149: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
151: PetscInt *work,*ia,*ja,nz,nslim_row,n,row,col,*ns_col,nslim_col;
152: PetscInt *tns,*tvc,nsz,i1,i2;
153: const PetscInt *j,*ai= a->i,*aj = a->j,*ns_row = a->inode.size;
156: nslim_row = a->inode.node_count;
157: n = A->cmap->n;
159: /* Create The column_inode for this matrix */
160: MatCreateColInode_Private(A,&nslim_col,&ns_col);
162: /* allocate space for reformated column_inode structure */
163: PetscMalloc1(nslim_col +1,&tns);
164: PetscMalloc1(n + 1,&tvc);
165: for (i1=0,tns[0]=0; i1<nslim_col; ++i1) tns[i1+1] = tns[i1] + ns_col[i1];
167: for (i1=0,col=0; i1<nslim_col; ++i1) {
168: nsz = ns_col[i1];
169: for (i2=0; i2<nsz; ++i2,++col) tvc[col] = i1;
170: }
171: /* allocate space for row pointers */
172: PetscMalloc1(nslim_row+1,&ia);
173: *iia = ia;
174: PetscMemzero(ia,(nslim_row+1)*sizeof(PetscInt));
175: PetscMalloc1(nslim_row+1,&work);
177: /* determine the number of columns in each row */
178: ia[0] = oshift;
179: for (i1=0,row=0; i1<nslim_row; row+=ns_row[i1],i1++) {
180: j = aj + ai[row] + ishift;
181: nz = ai[row+1] - ai[row];
182: if (!nz) continue; /* empty row */
183: col = *j++ + ishift;
184: i2 = tvc[col];
185: while (nz-- > 0) { /* off-diagonal elemets */
186: ia[i1+1]++;
187: i2++; /* Start col of next node */
188: while (nz > 0 && ((col = *j++ + ishift) < tns[i2])) nz--;
189: if (nz > 0) i2 = tvc[col];
190: }
191: }
193: /* shift ia[i] to point to next row */
194: for (i1=1; i1<nslim_row+1; i1++) {
195: row = ia[i1-1];
196: ia[i1] += row;
197: work[i1-1] = row - oshift;
198: }
200: /* allocate space for column pointers */
201: nz = ia[nslim_row] + (!ishift);
202: PetscMalloc1(nz,&ja);
203: *jja = ja;
205: /* loop over matrix putting into ja */
206: for (i1=0,row=0; i1<nslim_row; row+=ns_row[i1],i1++) {
207: j = aj + ai[row] + ishift;
208: nz = ai[row+1] - ai[row];
209: if (!nz) continue; /* empty row */
210: col = *j++ + ishift;
211: i2 = tvc[col];
212: while (nz-- > 0) {
213: ja[work[i1]++] = i2 + oshift;
214: ++i2;
215: while (nz > 0 && ((col = *j++ + ishift) < tns[i2])) nz--;
216: if (nz > 0) i2 = tvc[col];
217: }
218: }
219: PetscFree(ns_col);
220: PetscFree(work);
221: PetscFree(tns);
222: PetscFree(tvc);
223: return(0);
224: }
226: static PetscErrorCode MatGetRowIJ_SeqAIJ_Inode(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
227: {
228: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
232: *n = a->inode.node_count;
233: if (!ia) return(0);
234: if (!blockcompressed) {
235: MatGetRowIJ_SeqAIJ(A,oshift,symmetric,blockcompressed,n,ia,ja,done);;
236: } else if (symmetric) {
237: MatGetRowIJ_SeqAIJ_Inode_Symmetric(A,ia,ja,0,oshift);
238: } else {
239: MatGetRowIJ_SeqAIJ_Inode_Nonsymmetric(A,ia,ja,0,oshift);
240: }
241: return(0);
242: }
244: static PetscErrorCode MatRestoreRowIJ_SeqAIJ_Inode(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
245: {
249: if (!ia) return(0);
251: if (!blockcompressed) {
252: MatRestoreRowIJ_SeqAIJ(A,oshift,symmetric,blockcompressed,n,ia,ja,done);;
253: } else {
254: PetscFree(*ia);
255: PetscFree(*ja);
256: }
257: return(0);
258: }
260: /* ----------------------------------------------------------- */
262: static PetscErrorCode MatGetColumnIJ_SeqAIJ_Inode_Nonsymmetric(Mat A,const PetscInt *iia[],const PetscInt *jja[],PetscInt ishift,PetscInt oshift)
263: {
264: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
266: PetscInt *work,*ia,*ja,*j,nz,nslim_row, n,row,col,*ns_col,nslim_col;
267: PetscInt *tns,*tvc,*ns_row = a->inode.size,nsz,i1,i2,*ai= a->i,*aj = a->j;
270: nslim_row = a->inode.node_count;
271: n = A->cmap->n;
273: /* Create The column_inode for this matrix */
274: MatCreateColInode_Private(A,&nslim_col,&ns_col);
276: /* allocate space for reformated column_inode structure */
277: PetscMalloc1(nslim_col + 1,&tns);
278: PetscMalloc1(n + 1,&tvc);
279: for (i1=0,tns[0]=0; i1<nslim_col; ++i1) tns[i1+1] = tns[i1] + ns_col[i1];
281: for (i1=0,col=0; i1<nslim_col; ++i1) {
282: nsz = ns_col[i1];
283: for (i2=0; i2<nsz; ++i2,++col) tvc[col] = i1;
284: }
285: /* allocate space for column pointers */
286: PetscMalloc1(nslim_col+1,&ia);
287: *iia = ia;
288: PetscMemzero(ia,(nslim_col+1)*sizeof(PetscInt));
289: PetscMalloc1(nslim_col+1,&work);
291: /* determine the number of columns in each row */
292: ia[0] = oshift;
293: for (i1=0,row=0; i1<nslim_row; row+=ns_row[i1],i1++) {
294: j = aj + ai[row] + ishift;
295: col = *j++ + ishift;
296: i2 = tvc[col];
297: nz = ai[row+1] - ai[row];
298: while (nz-- > 0) { /* off-diagonal elemets */
299: /* ia[i1+1]++; */
300: ia[i2+1]++;
301: i2++;
302: while (nz > 0 && ((col = *j++ + ishift) < tns[i2])) nz--;
303: if (nz > 0) i2 = tvc[col];
304: }
305: }
307: /* shift ia[i] to point to next col */
308: for (i1=1; i1<nslim_col+1; i1++) {
309: col = ia[i1-1];
310: ia[i1] += col;
311: work[i1-1] = col - oshift;
312: }
314: /* allocate space for column pointers */
315: nz = ia[nslim_col] + (!ishift);
316: PetscMalloc1(nz,&ja);
317: *jja = ja;
319: /* loop over matrix putting into ja */
320: for (i1=0,row=0; i1<nslim_row; row+=ns_row[i1],i1++) {
321: j = aj + ai[row] + ishift;
322: col = *j++ + ishift;
323: i2 = tvc[col];
324: nz = ai[row+1] - ai[row];
325: while (nz-- > 0) {
326: /* ja[work[i1]++] = i2 + oshift; */
327: ja[work[i2]++] = i1 + oshift;
328: i2++;
329: while (nz > 0 && ((col = *j++ + ishift) < tns[i2])) nz--;
330: if (nz > 0) i2 = tvc[col];
331: }
332: }
333: PetscFree(ns_col);
334: PetscFree(work);
335: PetscFree(tns);
336: PetscFree(tvc);
337: return(0);
338: }
340: static PetscErrorCode MatGetColumnIJ_SeqAIJ_Inode(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
341: {
345: MatCreateColInode_Private(A,n,NULL);
346: if (!ia) return(0);
348: if (!blockcompressed) {
349: MatGetColumnIJ_SeqAIJ(A,oshift,symmetric,blockcompressed,n,ia,ja,done);;
350: } else if (symmetric) {
351: /* Since the indices are symmetric it does'nt matter */
352: MatGetRowIJ_SeqAIJ_Inode_Symmetric(A,ia,ja,0,oshift);
353: } else {
354: MatGetColumnIJ_SeqAIJ_Inode_Nonsymmetric(A,ia,ja,0,oshift);
355: }
356: return(0);
357: }
359: static PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Inode(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done)
360: {
364: if (!ia) return(0);
365: if (!blockcompressed) {
366: MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,blockcompressed,n,ia,ja,done);;
367: } else {
368: PetscFree(*ia);
369: PetscFree(*ja);
370: }
371: return(0);
372: }
374: /* ----------------------------------------------------------- */
376: static PetscErrorCode MatMult_SeqAIJ_Inode(Mat A,Vec xx,Vec yy)
377: {
378: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
379: PetscScalar sum1,sum2,sum3,sum4,sum5,tmp0,tmp1;
380: PetscScalar *y;
381: const PetscScalar *x;
382: const MatScalar *v1,*v2,*v3,*v4,*v5;
383: PetscErrorCode ierr;
384: PetscInt i1,i2,n,i,row,node_max,nsz,sz,nonzerorow=0;
385: const PetscInt *idx,*ns,*ii;
387: #if defined(PETSC_HAVE_PRAGMA_DISJOINT)
388: #pragma disjoint(*x,*y,*v1,*v2,*v3,*v4,*v5)
389: #endif
392: if (!a->inode.size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Missing Inode Structure");
393: node_max = a->inode.node_count;
394: ns = a->inode.size; /* Node Size array */
395: VecGetArrayRead(xx,&x);
396: VecGetArray(yy,&y);
397: idx = a->j;
398: v1 = a->a;
399: ii = a->i;
401: for (i = 0,row = 0; i< node_max; ++i) {
402: nsz = ns[i];
403: n = ii[1] - ii[0];
404: nonzerorow += (n>0)*nsz;
405: ii += nsz;
406: PetscPrefetchBlock(idx+nsz*n,n,0,PETSC_PREFETCH_HINT_NTA); /* Prefetch the indices for the block row after the current one */
407: PetscPrefetchBlock(v1+nsz*n,nsz*n,0,PETSC_PREFETCH_HINT_NTA); /* Prefetch the values for the block row after the current one */
408: sz = n; /* No of non zeros in this row */
409: /* Switch on the size of Node */
410: switch (nsz) { /* Each loop in 'case' is unrolled */
411: case 1:
412: sum1 = 0.;
414: for (n = 0; n< sz-1; n+=2) {
415: i1 = idx[0]; /* The instructions are ordered to */
416: i2 = idx[1]; /* make the compiler's job easy */
417: idx += 2;
418: tmp0 = x[i1];
419: tmp1 = x[i2];
420: sum1 += v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
421: }
423: if (n == sz-1) { /* Take care of the last nonzero */
424: tmp0 = x[*idx++];
425: sum1 += *v1++ *tmp0;
426: }
427: y[row++]=sum1;
428: break;
429: case 2:
430: sum1 = 0.;
431: sum2 = 0.;
432: v2 = v1 + n;
434: for (n = 0; n< sz-1; n+=2) {
435: i1 = idx[0];
436: i2 = idx[1];
437: idx += 2;
438: tmp0 = x[i1];
439: tmp1 = x[i2];
440: sum1 += v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
441: sum2 += v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
442: }
443: if (n == sz-1) {
444: tmp0 = x[*idx++];
445: sum1 += *v1++ * tmp0;
446: sum2 += *v2++ * tmp0;
447: }
448: y[row++]=sum1;
449: y[row++]=sum2;
450: v1 =v2; /* Since the next block to be processed starts there*/
451: idx +=sz;
452: break;
453: case 3:
454: sum1 = 0.;
455: sum2 = 0.;
456: sum3 = 0.;
457: v2 = v1 + n;
458: v3 = v2 + n;
460: for (n = 0; n< sz-1; n+=2) {
461: i1 = idx[0];
462: i2 = idx[1];
463: idx += 2;
464: tmp0 = x[i1];
465: tmp1 = x[i2];
466: sum1 += v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
467: sum2 += v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
468: sum3 += v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
469: }
470: if (n == sz-1) {
471: tmp0 = x[*idx++];
472: sum1 += *v1++ * tmp0;
473: sum2 += *v2++ * tmp0;
474: sum3 += *v3++ * tmp0;
475: }
476: y[row++]=sum1;
477: y[row++]=sum2;
478: y[row++]=sum3;
479: v1 =v3; /* Since the next block to be processed starts there*/
480: idx +=2*sz;
481: break;
482: case 4:
483: sum1 = 0.;
484: sum2 = 0.;
485: sum3 = 0.;
486: sum4 = 0.;
487: v2 = v1 + n;
488: v3 = v2 + n;
489: v4 = v3 + n;
491: for (n = 0; n< sz-1; n+=2) {
492: i1 = idx[0];
493: i2 = idx[1];
494: idx += 2;
495: tmp0 = x[i1];
496: tmp1 = x[i2];
497: sum1 += v1[0] * tmp0 + v1[1] *tmp1; v1 += 2;
498: sum2 += v2[0] * tmp0 + v2[1] *tmp1; v2 += 2;
499: sum3 += v3[0] * tmp0 + v3[1] *tmp1; v3 += 2;
500: sum4 += v4[0] * tmp0 + v4[1] *tmp1; v4 += 2;
501: }
502: if (n == sz-1) {
503: tmp0 = x[*idx++];
504: sum1 += *v1++ * tmp0;
505: sum2 += *v2++ * tmp0;
506: sum3 += *v3++ * tmp0;
507: sum4 += *v4++ * tmp0;
508: }
509: y[row++]=sum1;
510: y[row++]=sum2;
511: y[row++]=sum3;
512: y[row++]=sum4;
513: v1 =v4; /* Since the next block to be processed starts there*/
514: idx +=3*sz;
515: break;
516: case 5:
517: sum1 = 0.;
518: sum2 = 0.;
519: sum3 = 0.;
520: sum4 = 0.;
521: sum5 = 0.;
522: v2 = v1 + n;
523: v3 = v2 + n;
524: v4 = v3 + n;
525: v5 = v4 + n;
527: for (n = 0; n<sz-1; n+=2) {
528: i1 = idx[0];
529: i2 = idx[1];
530: idx += 2;
531: tmp0 = x[i1];
532: tmp1 = x[i2];
533: sum1 += v1[0] * tmp0 + v1[1] *tmp1; v1 += 2;
534: sum2 += v2[0] * tmp0 + v2[1] *tmp1; v2 += 2;
535: sum3 += v3[0] * tmp0 + v3[1] *tmp1; v3 += 2;
536: sum4 += v4[0] * tmp0 + v4[1] *tmp1; v4 += 2;
537: sum5 += v5[0] * tmp0 + v5[1] *tmp1; v5 += 2;
538: }
539: if (n == sz-1) {
540: tmp0 = x[*idx++];
541: sum1 += *v1++ * tmp0;
542: sum2 += *v2++ * tmp0;
543: sum3 += *v3++ * tmp0;
544: sum4 += *v4++ * tmp0;
545: sum5 += *v5++ * tmp0;
546: }
547: y[row++]=sum1;
548: y[row++]=sum2;
549: y[row++]=sum3;
550: y[row++]=sum4;
551: y[row++]=sum5;
552: v1 =v5; /* Since the next block to be processed starts there */
553: idx +=4*sz;
554: break;
555: default:
556: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Node size not yet supported");
557: }
558: }
559: VecRestoreArrayRead(xx,&x);
560: VecRestoreArray(yy,&y);
561: PetscLogFlops(2.0*a->nz - nonzerorow);
562: return(0);
563: }
564: /* ----------------------------------------------------------- */
565: /* Almost same code as the MatMult_SeqAIJ_Inode() */
566: static PetscErrorCode MatMultAdd_SeqAIJ_Inode(Mat A,Vec xx,Vec zz,Vec yy)
567: {
568: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
569: PetscScalar sum1,sum2,sum3,sum4,sum5,tmp0,tmp1;
570: const MatScalar *v1,*v2,*v3,*v4,*v5;
571: const PetscScalar *x;
572: PetscScalar *y,*z,*zt;
573: PetscErrorCode ierr;
574: PetscInt i1,i2,n,i,row,node_max,nsz,sz;
575: const PetscInt *idx,*ns,*ii;
578: if (!a->inode.size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Missing Inode Structure");
579: node_max = a->inode.node_count;
580: ns = a->inode.size; /* Node Size array */
582: VecGetArrayRead(xx,&x);
583: VecGetArrayPair(zz,yy,&z,&y);
584: zt = z;
586: idx = a->j;
587: v1 = a->a;
588: ii = a->i;
590: for (i = 0,row = 0; i< node_max; ++i) {
591: nsz = ns[i];
592: n = ii[1] - ii[0];
593: ii += nsz;
594: sz = n; /* No of non zeros in this row */
595: /* Switch on the size of Node */
596: switch (nsz) { /* Each loop in 'case' is unrolled */
597: case 1:
598: sum1 = *zt++;
600: for (n = 0; n< sz-1; n+=2) {
601: i1 = idx[0]; /* The instructions are ordered to */
602: i2 = idx[1]; /* make the compiler's job easy */
603: idx += 2;
604: tmp0 = x[i1];
605: tmp1 = x[i2];
606: sum1 += v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
607: }
609: if (n == sz-1) { /* Take care of the last nonzero */
610: tmp0 = x[*idx++];
611: sum1 += *v1++ * tmp0;
612: }
613: y[row++]=sum1;
614: break;
615: case 2:
616: sum1 = *zt++;
617: sum2 = *zt++;
618: v2 = v1 + n;
620: for (n = 0; n< sz-1; n+=2) {
621: i1 = idx[0];
622: i2 = idx[1];
623: idx += 2;
624: tmp0 = x[i1];
625: tmp1 = x[i2];
626: sum1 += v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
627: sum2 += v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
628: }
629: if (n == sz-1) {
630: tmp0 = x[*idx++];
631: sum1 += *v1++ * tmp0;
632: sum2 += *v2++ * tmp0;
633: }
634: y[row++]=sum1;
635: y[row++]=sum2;
636: v1 =v2; /* Since the next block to be processed starts there*/
637: idx +=sz;
638: break;
639: case 3:
640: sum1 = *zt++;
641: sum2 = *zt++;
642: sum3 = *zt++;
643: v2 = v1 + n;
644: v3 = v2 + n;
646: for (n = 0; n< sz-1; n+=2) {
647: i1 = idx[0];
648: i2 = idx[1];
649: idx += 2;
650: tmp0 = x[i1];
651: tmp1 = x[i2];
652: sum1 += v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
653: sum2 += v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
654: sum3 += v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
655: }
656: if (n == sz-1) {
657: tmp0 = x[*idx++];
658: sum1 += *v1++ * tmp0;
659: sum2 += *v2++ * tmp0;
660: sum3 += *v3++ * tmp0;
661: }
662: y[row++]=sum1;
663: y[row++]=sum2;
664: y[row++]=sum3;
665: v1 =v3; /* Since the next block to be processed starts there*/
666: idx +=2*sz;
667: break;
668: case 4:
669: sum1 = *zt++;
670: sum2 = *zt++;
671: sum3 = *zt++;
672: sum4 = *zt++;
673: v2 = v1 + n;
674: v3 = v2 + n;
675: v4 = v3 + n;
677: for (n = 0; n< sz-1; n+=2) {
678: i1 = idx[0];
679: i2 = idx[1];
680: idx += 2;
681: tmp0 = x[i1];
682: tmp1 = x[i2];
683: sum1 += v1[0] * tmp0 + v1[1] *tmp1; v1 += 2;
684: sum2 += v2[0] * tmp0 + v2[1] *tmp1; v2 += 2;
685: sum3 += v3[0] * tmp0 + v3[1] *tmp1; v3 += 2;
686: sum4 += v4[0] * tmp0 + v4[1] *tmp1; v4 += 2;
687: }
688: if (n == sz-1) {
689: tmp0 = x[*idx++];
690: sum1 += *v1++ * tmp0;
691: sum2 += *v2++ * tmp0;
692: sum3 += *v3++ * tmp0;
693: sum4 += *v4++ * tmp0;
694: }
695: y[row++]=sum1;
696: y[row++]=sum2;
697: y[row++]=sum3;
698: y[row++]=sum4;
699: v1 =v4; /* Since the next block to be processed starts there*/
700: idx +=3*sz;
701: break;
702: case 5:
703: sum1 = *zt++;
704: sum2 = *zt++;
705: sum3 = *zt++;
706: sum4 = *zt++;
707: sum5 = *zt++;
708: v2 = v1 + n;
709: v3 = v2 + n;
710: v4 = v3 + n;
711: v5 = v4 + n;
713: for (n = 0; n<sz-1; n+=2) {
714: i1 = idx[0];
715: i2 = idx[1];
716: idx += 2;
717: tmp0 = x[i1];
718: tmp1 = x[i2];
719: sum1 += v1[0] * tmp0 + v1[1] *tmp1; v1 += 2;
720: sum2 += v2[0] * tmp0 + v2[1] *tmp1; v2 += 2;
721: sum3 += v3[0] * tmp0 + v3[1] *tmp1; v3 += 2;
722: sum4 += v4[0] * tmp0 + v4[1] *tmp1; v4 += 2;
723: sum5 += v5[0] * tmp0 + v5[1] *tmp1; v5 += 2;
724: }
725: if (n == sz-1) {
726: tmp0 = x[*idx++];
727: sum1 += *v1++ * tmp0;
728: sum2 += *v2++ * tmp0;
729: sum3 += *v3++ * tmp0;
730: sum4 += *v4++ * tmp0;
731: sum5 += *v5++ * tmp0;
732: }
733: y[row++]=sum1;
734: y[row++]=sum2;
735: y[row++]=sum3;
736: y[row++]=sum4;
737: y[row++]=sum5;
738: v1 =v5; /* Since the next block to be processed starts there */
739: idx +=4*sz;
740: break;
741: default:
742: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Node size not yet supported");
743: }
744: }
745: VecRestoreArrayRead(xx,&x);
746: VecRestoreArrayPair(zz,yy,&z,&y);
747: PetscLogFlops(2.0*a->nz);
748: return(0);
749: }
751: /* ----------------------------------------------------------- */
752: PetscErrorCode MatSolve_SeqAIJ_Inode_inplace(Mat A,Vec bb,Vec xx)
753: {
754: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
755: IS iscol = a->col,isrow = a->row;
756: PetscErrorCode ierr;
757: const PetscInt *r,*c,*rout,*cout;
758: PetscInt i,j,n = A->rmap->n,nz;
759: PetscInt node_max,*ns,row,nsz,aii,i0,i1;
760: const PetscInt *ai = a->i,*a_j = a->j,*vi,*ad,*aj;
761: PetscScalar *x,*tmp,*tmps,tmp0,tmp1;
762: PetscScalar sum1,sum2,sum3,sum4,sum5;
763: const MatScalar *v1,*v2,*v3,*v4,*v5,*a_a = a->a,*aa;
764: const PetscScalar *b;
767: if (!a->inode.size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Missing Inode Structure");
768: node_max = a->inode.node_count;
769: ns = a->inode.size; /* Node Size array */
771: VecGetArrayRead(bb,&b);
772: VecGetArray(xx,&x);
773: tmp = a->solve_work;
775: ISGetIndices(isrow,&rout); r = rout;
776: ISGetIndices(iscol,&cout); c = cout + (n-1);
778: /* forward solve the lower triangular */
779: tmps = tmp;
780: aa = a_a;
781: aj = a_j;
782: ad = a->diag;
784: for (i = 0,row = 0; i< node_max; ++i) {
785: nsz = ns[i];
786: aii = ai[row];
787: v1 = aa + aii;
788: vi = aj + aii;
789: nz = ad[row]- aii;
790: if (i < node_max-1) {
791: /* Prefetch the block after the current one, the prefetch itself can't cause a memory error,
792: * but our indexing to determine it's size could. */
793: PetscPrefetchBlock(aj+ai[row+nsz],ad[row+nsz]-ai[row+nsz],0,PETSC_PREFETCH_HINT_NTA); /* indices */
794: /* In my tests, it seems to be better to fetch entire rows instead of just the lower-triangular part */
795: PetscPrefetchBlock(aa+ai[row+nsz],ad[row+nsz+ns[i+1]-1]-ai[row+nsz],0,PETSC_PREFETCH_HINT_NTA);
796: /* for (j=0; j<ns[i+1]; j++) PetscPrefetchBlock(aa+ai[row+nsz+j],ad[row+nsz+j]-ai[row+nsz+j],0,0); */
797: }
799: switch (nsz) { /* Each loop in 'case' is unrolled */
800: case 1:
801: sum1 = b[*r++];
802: for (j=0; j<nz-1; j+=2) {
803: i0 = vi[0];
804: i1 = vi[1];
805: vi +=2;
806: tmp0 = tmps[i0];
807: tmp1 = tmps[i1];
808: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
809: }
810: if (j == nz-1) {
811: tmp0 = tmps[*vi++];
812: sum1 -= *v1++ *tmp0;
813: }
814: tmp[row++]=sum1;
815: break;
816: case 2:
817: sum1 = b[*r++];
818: sum2 = b[*r++];
819: v2 = aa + ai[row+1];
821: for (j=0; j<nz-1; j+=2) {
822: i0 = vi[0];
823: i1 = vi[1];
824: vi +=2;
825: tmp0 = tmps[i0];
826: tmp1 = tmps[i1];
827: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
828: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
829: }
830: if (j == nz-1) {
831: tmp0 = tmps[*vi++];
832: sum1 -= *v1++ *tmp0;
833: sum2 -= *v2++ *tmp0;
834: }
835: sum2 -= *v2++ *sum1;
836: tmp[row++]=sum1;
837: tmp[row++]=sum2;
838: break;
839: case 3:
840: sum1 = b[*r++];
841: sum2 = b[*r++];
842: sum3 = b[*r++];
843: v2 = aa + ai[row+1];
844: v3 = aa + ai[row+2];
846: for (j=0; j<nz-1; j+=2) {
847: i0 = vi[0];
848: i1 = vi[1];
849: vi +=2;
850: tmp0 = tmps[i0];
851: tmp1 = tmps[i1];
852: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
853: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
854: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
855: }
856: if (j == nz-1) {
857: tmp0 = tmps[*vi++];
858: sum1 -= *v1++ *tmp0;
859: sum2 -= *v2++ *tmp0;
860: sum3 -= *v3++ *tmp0;
861: }
862: sum2 -= *v2++ * sum1;
863: sum3 -= *v3++ * sum1;
864: sum3 -= *v3++ * sum2;
866: tmp[row++]=sum1;
867: tmp[row++]=sum2;
868: tmp[row++]=sum3;
869: break;
871: case 4:
872: sum1 = b[*r++];
873: sum2 = b[*r++];
874: sum3 = b[*r++];
875: sum4 = b[*r++];
876: v2 = aa + ai[row+1];
877: v3 = aa + ai[row+2];
878: v4 = aa + ai[row+3];
880: for (j=0; j<nz-1; j+=2) {
881: i0 = vi[0];
882: i1 = vi[1];
883: vi +=2;
884: tmp0 = tmps[i0];
885: tmp1 = tmps[i1];
886: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
887: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
888: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
889: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
890: }
891: if (j == nz-1) {
892: tmp0 = tmps[*vi++];
893: sum1 -= *v1++ *tmp0;
894: sum2 -= *v2++ *tmp0;
895: sum3 -= *v3++ *tmp0;
896: sum4 -= *v4++ *tmp0;
897: }
898: sum2 -= *v2++ * sum1;
899: sum3 -= *v3++ * sum1;
900: sum4 -= *v4++ * sum1;
901: sum3 -= *v3++ * sum2;
902: sum4 -= *v4++ * sum2;
903: sum4 -= *v4++ * sum3;
905: tmp[row++]=sum1;
906: tmp[row++]=sum2;
907: tmp[row++]=sum3;
908: tmp[row++]=sum4;
909: break;
910: case 5:
911: sum1 = b[*r++];
912: sum2 = b[*r++];
913: sum3 = b[*r++];
914: sum4 = b[*r++];
915: sum5 = b[*r++];
916: v2 = aa + ai[row+1];
917: v3 = aa + ai[row+2];
918: v4 = aa + ai[row+3];
919: v5 = aa + ai[row+4];
921: for (j=0; j<nz-1; j+=2) {
922: i0 = vi[0];
923: i1 = vi[1];
924: vi +=2;
925: tmp0 = tmps[i0];
926: tmp1 = tmps[i1];
927: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
928: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
929: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
930: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
931: sum5 -= v5[0] * tmp0 + v5[1] * tmp1; v5 += 2;
932: }
933: if (j == nz-1) {
934: tmp0 = tmps[*vi++];
935: sum1 -= *v1++ *tmp0;
936: sum2 -= *v2++ *tmp0;
937: sum3 -= *v3++ *tmp0;
938: sum4 -= *v4++ *tmp0;
939: sum5 -= *v5++ *tmp0;
940: }
942: sum2 -= *v2++ * sum1;
943: sum3 -= *v3++ * sum1;
944: sum4 -= *v4++ * sum1;
945: sum5 -= *v5++ * sum1;
946: sum3 -= *v3++ * sum2;
947: sum4 -= *v4++ * sum2;
948: sum5 -= *v5++ * sum2;
949: sum4 -= *v4++ * sum3;
950: sum5 -= *v5++ * sum3;
951: sum5 -= *v5++ * sum4;
953: tmp[row++]=sum1;
954: tmp[row++]=sum2;
955: tmp[row++]=sum3;
956: tmp[row++]=sum4;
957: tmp[row++]=sum5;
958: break;
959: default:
960: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Node size not yet supported \n");
961: }
962: }
963: /* backward solve the upper triangular */
964: for (i=node_max -1,row = n-1; i>=0; i--) {
965: nsz = ns[i];
966: aii = ai[row+1] -1;
967: v1 = aa + aii;
968: vi = aj + aii;
969: nz = aii- ad[row];
970: switch (nsz) { /* Each loop in 'case' is unrolled */
971: case 1:
972: sum1 = tmp[row];
974: for (j=nz; j>1; j-=2) {
975: vi -=2;
976: i0 = vi[2];
977: i1 = vi[1];
978: tmp0 = tmps[i0];
979: tmp1 = tmps[i1];
980: v1 -= 2;
981: sum1 -= v1[2] * tmp0 + v1[1] * tmp1;
982: }
983: if (j==1) {
984: tmp0 = tmps[*vi--];
985: sum1 -= *v1-- * tmp0;
986: }
987: x[*c--] = tmp[row] = sum1*a_a[ad[row]]; row--;
988: break;
989: case 2:
990: sum1 = tmp[row];
991: sum2 = tmp[row -1];
992: v2 = aa + ai[row]-1;
993: for (j=nz; j>1; j-=2) {
994: vi -=2;
995: i0 = vi[2];
996: i1 = vi[1];
997: tmp0 = tmps[i0];
998: tmp1 = tmps[i1];
999: v1 -= 2;
1000: v2 -= 2;
1001: sum1 -= v1[2] * tmp0 + v1[1] * tmp1;
1002: sum2 -= v2[2] * tmp0 + v2[1] * tmp1;
1003: }
1004: if (j==1) {
1005: tmp0 = tmps[*vi--];
1006: sum1 -= *v1-- * tmp0;
1007: sum2 -= *v2-- * tmp0;
1008: }
1010: tmp0 = x[*c--] = tmp[row] = sum1*a_a[ad[row]]; row--;
1011: sum2 -= *v2-- * tmp0;
1012: x[*c--] = tmp[row] = sum2*a_a[ad[row]]; row--;
1013: break;
1014: case 3:
1015: sum1 = tmp[row];
1016: sum2 = tmp[row -1];
1017: sum3 = tmp[row -2];
1018: v2 = aa + ai[row]-1;
1019: v3 = aa + ai[row -1]-1;
1020: for (j=nz; j>1; j-=2) {
1021: vi -=2;
1022: i0 = vi[2];
1023: i1 = vi[1];
1024: tmp0 = tmps[i0];
1025: tmp1 = tmps[i1];
1026: v1 -= 2;
1027: v2 -= 2;
1028: v3 -= 2;
1029: sum1 -= v1[2] * tmp0 + v1[1] * tmp1;
1030: sum2 -= v2[2] * tmp0 + v2[1] * tmp1;
1031: sum3 -= v3[2] * tmp0 + v3[1] * tmp1;
1032: }
1033: if (j==1) {
1034: tmp0 = tmps[*vi--];
1035: sum1 -= *v1-- * tmp0;
1036: sum2 -= *v2-- * tmp0;
1037: sum3 -= *v3-- * tmp0;
1038: }
1039: tmp0 = x[*c--] = tmp[row] = sum1*a_a[ad[row]]; row--;
1040: sum2 -= *v2-- * tmp0;
1041: sum3 -= *v3-- * tmp0;
1042: tmp0 = x[*c--] = tmp[row] = sum2*a_a[ad[row]]; row--;
1043: sum3 -= *v3-- * tmp0;
1044: x[*c--] = tmp[row] = sum3*a_a[ad[row]]; row--;
1046: break;
1047: case 4:
1048: sum1 = tmp[row];
1049: sum2 = tmp[row -1];
1050: sum3 = tmp[row -2];
1051: sum4 = tmp[row -3];
1052: v2 = aa + ai[row]-1;
1053: v3 = aa + ai[row -1]-1;
1054: v4 = aa + ai[row -2]-1;
1056: for (j=nz; j>1; j-=2) {
1057: vi -=2;
1058: i0 = vi[2];
1059: i1 = vi[1];
1060: tmp0 = tmps[i0];
1061: tmp1 = tmps[i1];
1062: v1 -= 2;
1063: v2 -= 2;
1064: v3 -= 2;
1065: v4 -= 2;
1066: sum1 -= v1[2] * tmp0 + v1[1] * tmp1;
1067: sum2 -= v2[2] * tmp0 + v2[1] * tmp1;
1068: sum3 -= v3[2] * tmp0 + v3[1] * tmp1;
1069: sum4 -= v4[2] * tmp0 + v4[1] * tmp1;
1070: }
1071: if (j==1) {
1072: tmp0 = tmps[*vi--];
1073: sum1 -= *v1-- * tmp0;
1074: sum2 -= *v2-- * tmp0;
1075: sum3 -= *v3-- * tmp0;
1076: sum4 -= *v4-- * tmp0;
1077: }
1079: tmp0 = x[*c--] = tmp[row] = sum1*a_a[ad[row]]; row--;
1080: sum2 -= *v2-- * tmp0;
1081: sum3 -= *v3-- * tmp0;
1082: sum4 -= *v4-- * tmp0;
1083: tmp0 = x[*c--] = tmp[row] = sum2*a_a[ad[row]]; row--;
1084: sum3 -= *v3-- * tmp0;
1085: sum4 -= *v4-- * tmp0;
1086: tmp0 = x[*c--] = tmp[row] = sum3*a_a[ad[row]]; row--;
1087: sum4 -= *v4-- * tmp0;
1088: x[*c--] = tmp[row] = sum4*a_a[ad[row]]; row--;
1089: break;
1090: case 5:
1091: sum1 = tmp[row];
1092: sum2 = tmp[row -1];
1093: sum3 = tmp[row -2];
1094: sum4 = tmp[row -3];
1095: sum5 = tmp[row -4];
1096: v2 = aa + ai[row]-1;
1097: v3 = aa + ai[row -1]-1;
1098: v4 = aa + ai[row -2]-1;
1099: v5 = aa + ai[row -3]-1;
1100: for (j=nz; j>1; j-=2) {
1101: vi -= 2;
1102: i0 = vi[2];
1103: i1 = vi[1];
1104: tmp0 = tmps[i0];
1105: tmp1 = tmps[i1];
1106: v1 -= 2;
1107: v2 -= 2;
1108: v3 -= 2;
1109: v4 -= 2;
1110: v5 -= 2;
1111: sum1 -= v1[2] * tmp0 + v1[1] * tmp1;
1112: sum2 -= v2[2] * tmp0 + v2[1] * tmp1;
1113: sum3 -= v3[2] * tmp0 + v3[1] * tmp1;
1114: sum4 -= v4[2] * tmp0 + v4[1] * tmp1;
1115: sum5 -= v5[2] * tmp0 + v5[1] * tmp1;
1116: }
1117: if (j==1) {
1118: tmp0 = tmps[*vi--];
1119: sum1 -= *v1-- * tmp0;
1120: sum2 -= *v2-- * tmp0;
1121: sum3 -= *v3-- * tmp0;
1122: sum4 -= *v4-- * tmp0;
1123: sum5 -= *v5-- * tmp0;
1124: }
1126: tmp0 = x[*c--] = tmp[row] = sum1*a_a[ad[row]]; row--;
1127: sum2 -= *v2-- * tmp0;
1128: sum3 -= *v3-- * tmp0;
1129: sum4 -= *v4-- * tmp0;
1130: sum5 -= *v5-- * tmp0;
1131: tmp0 = x[*c--] = tmp[row] = sum2*a_a[ad[row]]; row--;
1132: sum3 -= *v3-- * tmp0;
1133: sum4 -= *v4-- * tmp0;
1134: sum5 -= *v5-- * tmp0;
1135: tmp0 = x[*c--] = tmp[row] = sum3*a_a[ad[row]]; row--;
1136: sum4 -= *v4-- * tmp0;
1137: sum5 -= *v5-- * tmp0;
1138: tmp0 = x[*c--] = tmp[row] = sum4*a_a[ad[row]]; row--;
1139: sum5 -= *v5-- * tmp0;
1140: x[*c--] = tmp[row] = sum5*a_a[ad[row]]; row--;
1141: break;
1142: default:
1143: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Node size not yet supported \n");
1144: }
1145: }
1146: ISRestoreIndices(isrow,&rout);
1147: ISRestoreIndices(iscol,&cout);
1148: VecRestoreArrayRead(bb,&b);
1149: VecRestoreArray(xx,&x);
1150: PetscLogFlops(2.0*a->nz - A->cmap->n);
1151: return(0);
1152: }
1154: PetscErrorCode MatLUFactorNumeric_SeqAIJ_Inode(Mat B,Mat A,const MatFactorInfo *info)
1155: {
1156: Mat C =B;
1157: Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)C->data;
1158: IS isrow = b->row,isicol = b->icol;
1159: PetscErrorCode ierr;
1160: const PetscInt *r,*ic,*ics;
1161: const PetscInt n=A->rmap->n,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bdiag=b->diag;
1162: PetscInt i,j,k,nz,nzL,row,*pj;
1163: const PetscInt *ajtmp,*bjtmp;
1164: MatScalar *pc,*pc1,*pc2,*pc3,*pc4,mul1,mul2,mul3,mul4,*pv,*rtmp1,*rtmp2,*rtmp3,*rtmp4;
1165: const MatScalar *aa=a->a,*v,*v1,*v2,*v3,*v4;
1166: FactorShiftCtx sctx;
1167: const PetscInt *ddiag;
1168: PetscReal rs;
1169: MatScalar d;
1170: PetscInt inod,nodesz,node_max,col;
1171: const PetscInt *ns;
1172: PetscInt *tmp_vec1,*tmp_vec2,*nsmap;
1175: /* MatPivotSetUp(): initialize shift context sctx */
1176: PetscMemzero(&sctx,sizeof(FactorShiftCtx));
1178: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1179: ddiag = a->diag;
1180: sctx.shift_top = info->zeropivot;
1181: for (i=0; i<n; i++) {
1182: /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1183: d = (aa)[ddiag[i]];
1184: rs = -PetscAbsScalar(d) - PetscRealPart(d);
1185: v = aa+ai[i];
1186: nz = ai[i+1] - ai[i];
1187: for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]);
1188: if (rs>sctx.shift_top) sctx.shift_top = rs;
1189: }
1190: sctx.shift_top *= 1.1;
1191: sctx.nshift_max = 5;
1192: sctx.shift_lo = 0.;
1193: sctx.shift_hi = 1.;
1194: }
1196: ISGetIndices(isrow,&r);
1197: ISGetIndices(isicol,&ic);
1199: PetscCalloc4(n,&rtmp1,n,&rtmp2,n,&rtmp3,n,&rtmp4);
1200: ics = ic;
1202: node_max = a->inode.node_count;
1203: ns = a->inode.size;
1204: if (!ns) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Matrix without inode information");
1206: /* If max inode size > 4, split it into two inodes.*/
1207: /* also map the inode sizes according to the ordering */
1208: PetscMalloc1(n+1,&tmp_vec1);
1209: for (i=0,j=0; i<node_max; ++i,++j) {
1210: if (ns[i] > 4) {
1211: tmp_vec1[j] = 4;
1212: ++j;
1213: tmp_vec1[j] = ns[i] - tmp_vec1[j-1];
1214: } else {
1215: tmp_vec1[j] = ns[i];
1216: }
1217: }
1218: /* Use the correct node_max */
1219: node_max = j;
1221: /* Now reorder the inode info based on mat re-ordering info */
1222: /* First create a row -> inode_size_array_index map */
1223: PetscMalloc1(n+1,&nsmap);
1224: PetscMalloc1(node_max+1,&tmp_vec2);
1225: for (i=0,row=0; i<node_max; i++) {
1226: nodesz = tmp_vec1[i];
1227: for (j=0; j<nodesz; j++,row++) {
1228: nsmap[row] = i;
1229: }
1230: }
1231: /* Using nsmap, create a reordered ns structure */
1232: for (i=0,j=0; i< node_max; i++) {
1233: nodesz = tmp_vec1[nsmap[r[j]]]; /* here the reordered row_no is in r[] */
1234: tmp_vec2[i] = nodesz;
1235: j += nodesz;
1236: }
1237: PetscFree(nsmap);
1238: PetscFree(tmp_vec1);
1240: /* Now use the correct ns */
1241: ns = tmp_vec2;
1243: do {
1244: sctx.newshift = PETSC_FALSE;
1245: /* Now loop over each block-row, and do the factorization */
1246: for (inod=0,i=0; inod<node_max; inod++) { /* i: row index; inod: inode index */
1247: nodesz = ns[inod];
1249: switch (nodesz) {
1250: case 1:
1251: /*----------*/
1252: /* zero rtmp1 */
1253: /* L part */
1254: nz = bi[i+1] - bi[i];
1255: bjtmp = bj + bi[i];
1256: for (j=0; j<nz; j++) rtmp1[bjtmp[j]] = 0.0;
1258: /* U part */
1259: nz = bdiag[i]-bdiag[i+1];
1260: bjtmp = bj + bdiag[i+1]+1;
1261: for (j=0; j<nz; j++) rtmp1[bjtmp[j]] = 0.0;
1263: /* load in initial (unfactored row) */
1264: nz = ai[r[i]+1] - ai[r[i]];
1265: ajtmp = aj + ai[r[i]];
1266: v = aa + ai[r[i]];
1267: for (j=0; j<nz; j++) rtmp1[ics[ajtmp[j]]] = v[j];
1269: /* ZeropivotApply() */
1270: rtmp1[i] += sctx.shift_amount; /* shift the diagonal of the matrix */
1272: /* elimination */
1273: bjtmp = bj + bi[i];
1274: row = *bjtmp++;
1275: nzL = bi[i+1] - bi[i];
1276: for (k=0; k < nzL; k++) {
1277: pc = rtmp1 + row;
1278: if (*pc != 0.0) {
1279: pv = b->a + bdiag[row];
1280: mul1 = *pc * (*pv);
1281: *pc = mul1;
1282: pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
1283: pv = b->a + bdiag[row+1]+1;
1284: nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */
1285: for (j=0; j<nz; j++) rtmp1[pj[j]] -= mul1 * pv[j];
1286: PetscLogFlops(1+2*nz);
1287: }
1288: row = *bjtmp++;
1289: }
1291: /* finished row so stick it into b->a */
1292: rs = 0.0;
1293: /* L part */
1294: pv = b->a + bi[i];
1295: pj = b->j + bi[i];
1296: nz = bi[i+1] - bi[i];
1297: for (j=0; j<nz; j++) {
1298: pv[j] = rtmp1[pj[j]]; rs += PetscAbsScalar(pv[j]);
1299: }
1301: /* U part */
1302: pv = b->a + bdiag[i+1]+1;
1303: pj = b->j + bdiag[i+1]+1;
1304: nz = bdiag[i] - bdiag[i+1]-1;
1305: for (j=0; j<nz; j++) {
1306: pv[j] = rtmp1[pj[j]]; rs += PetscAbsScalar(pv[j]);
1307: }
1309: /* Check zero pivot */
1310: sctx.rs = rs;
1311: sctx.pv = rtmp1[i];
1312: MatPivotCheck(B,A,info,&sctx,i);
1313: if (sctx.newshift) break;
1315: /* Mark diagonal and invert diagonal for simplier triangular solves */
1316: pv = b->a + bdiag[i];
1317: *pv = 1.0/sctx.pv; /* sctx.pv = rtmp1[i]+shiftamount if shifttype==MAT_SHIFT_INBLOCKS */
1318: break;
1320: case 2:
1321: /*----------*/
1322: /* zero rtmp1 and rtmp2 */
1323: /* L part */
1324: nz = bi[i+1] - bi[i];
1325: bjtmp = bj + bi[i];
1326: for (j=0; j<nz; j++) {
1327: col = bjtmp[j];
1328: rtmp1[col] = 0.0; rtmp2[col] = 0.0;
1329: }
1331: /* U part */
1332: nz = bdiag[i]-bdiag[i+1];
1333: bjtmp = bj + bdiag[i+1]+1;
1334: for (j=0; j<nz; j++) {
1335: col = bjtmp[j];
1336: rtmp1[col] = 0.0; rtmp2[col] = 0.0;
1337: }
1339: /* load in initial (unfactored row) */
1340: nz = ai[r[i]+1] - ai[r[i]];
1341: ajtmp = aj + ai[r[i]];
1342: v1 = aa + ai[r[i]]; v2 = aa + ai[r[i]+1];
1343: for (j=0; j<nz; j++) {
1344: col = ics[ajtmp[j]];
1345: rtmp1[col] = v1[j]; rtmp2[col] = v2[j];
1346: }
1347: /* ZeropivotApply(): shift the diagonal of the matrix */
1348: rtmp1[i] += sctx.shift_amount; rtmp2[i+1] += sctx.shift_amount;
1350: /* elimination */
1351: bjtmp = bj + bi[i];
1352: row = *bjtmp++; /* pivot row */
1353: nzL = bi[i+1] - bi[i];
1354: for (k=0; k < nzL; k++) {
1355: pc1 = rtmp1 + row;
1356: pc2 = rtmp2 + row;
1357: if (*pc1 != 0.0 || *pc2 != 0.0) {
1358: pv = b->a + bdiag[row];
1359: mul1 = *pc1*(*pv); mul2 = *pc2*(*pv);
1360: *pc1 = mul1; *pc2 = mul2;
1362: pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
1363: pv = b->a + bdiag[row+1]+1;
1364: nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */
1365: for (j=0; j<nz; j++) {
1366: col = pj[j];
1367: rtmp1[col] -= mul1 * pv[j];
1368: rtmp2[col] -= mul2 * pv[j];
1369: }
1370: PetscLogFlops(2+4*nz);
1371: }
1372: row = *bjtmp++;
1373: }
1375: /* finished row i; check zero pivot, then stick row i into b->a */
1376: rs = 0.0;
1377: /* L part */
1378: pc1 = b->a + bi[i];
1379: pj = b->j + bi[i];
1380: nz = bi[i+1] - bi[i];
1381: for (j=0; j<nz; j++) {
1382: col = pj[j];
1383: pc1[j] = rtmp1[col]; rs += PetscAbsScalar(pc1[j]);
1384: }
1385: /* U part */
1386: pc1 = b->a + bdiag[i+1]+1;
1387: pj = b->j + bdiag[i+1]+1;
1388: nz = bdiag[i] - bdiag[i+1] - 1; /* exclude diagonal */
1389: for (j=0; j<nz; j++) {
1390: col = pj[j];
1391: pc1[j] = rtmp1[col]; rs += PetscAbsScalar(pc1[j]);
1392: }
1394: sctx.rs = rs;
1395: sctx.pv = rtmp1[i];
1396: MatPivotCheck(B,A,info,&sctx,i);
1397: if (sctx.newshift) break;
1398: pc1 = b->a + bdiag[i]; /* Mark diagonal */
1399: *pc1 = 1.0/sctx.pv;
1401: /* Now take care of diagonal 2x2 block. */
1402: pc2 = rtmp2 + i;
1403: if (*pc2 != 0.0) {
1404: mul1 = (*pc2)*(*pc1); /* *pc1=diag[i] is inverted! */
1405: *pc2 = mul1; /* insert L entry */
1406: pj = b->j + bdiag[i+1]+1; /* beginning of U(i,:) */
1407: nz = bdiag[i]-bdiag[i+1]-1; /* num of entries in U(i,:) excluding diag */
1408: for (j=0; j<nz; j++) {
1409: col = pj[j]; rtmp2[col] -= mul1 * rtmp1[col];
1410: }
1411: PetscLogFlops(1+2*nz);
1412: }
1414: /* finished row i+1; check zero pivot, then stick row i+1 into b->a */
1415: rs = 0.0;
1416: /* L part */
1417: pc2 = b->a + bi[i+1];
1418: pj = b->j + bi[i+1];
1419: nz = bi[i+2] - bi[i+1];
1420: for (j=0; j<nz; j++) {
1421: col = pj[j];
1422: pc2[j] = rtmp2[col]; rs += PetscAbsScalar(pc2[j]);
1423: }
1424: /* U part */
1425: pc2 = b->a + bdiag[i+2]+1;
1426: pj = b->j + bdiag[i+2]+1;
1427: nz = bdiag[i+1] - bdiag[i+2] - 1; /* exclude diagonal */
1428: for (j=0; j<nz; j++) {
1429: col = pj[j];
1430: pc2[j] = rtmp2[col]; rs += PetscAbsScalar(pc2[j]);
1431: }
1433: sctx.rs = rs;
1434: sctx.pv = rtmp2[i+1];
1435: MatPivotCheck(B,A,info,&sctx,i+1);
1436: if (sctx.newshift) break;
1437: pc2 = b->a + bdiag[i+1];
1438: *pc2 = 1.0/sctx.pv;
1439: break;
1441: case 3:
1442: /*----------*/
1443: /* zero rtmp */
1444: /* L part */
1445: nz = bi[i+1] - bi[i];
1446: bjtmp = bj + bi[i];
1447: for (j=0; j<nz; j++) {
1448: col = bjtmp[j];
1449: rtmp1[col] = 0.0; rtmp2[col] = 0.0; rtmp3[col] = 0.0;
1450: }
1452: /* U part */
1453: nz = bdiag[i]-bdiag[i+1];
1454: bjtmp = bj + bdiag[i+1]+1;
1455: for (j=0; j<nz; j++) {
1456: col = bjtmp[j];
1457: rtmp1[col] = 0.0; rtmp2[col] = 0.0; rtmp3[col] = 0.0;
1458: }
1460: /* load in initial (unfactored row) */
1461: nz = ai[r[i]+1] - ai[r[i]];
1462: ajtmp = aj + ai[r[i]];
1463: v1 = aa + ai[r[i]]; v2 = aa + ai[r[i]+1]; v3 = aa + ai[r[i]+2];
1464: for (j=0; j<nz; j++) {
1465: col = ics[ajtmp[j]];
1466: rtmp1[col] = v1[j]; rtmp2[col] = v2[j]; rtmp3[col] = v3[j];
1467: }
1468: /* ZeropivotApply(): shift the diagonal of the matrix */
1469: rtmp1[i] += sctx.shift_amount; rtmp2[i+1] += sctx.shift_amount; rtmp3[i+2] += sctx.shift_amount;
1471: /* elimination */
1472: bjtmp = bj + bi[i];
1473: row = *bjtmp++; /* pivot row */
1474: nzL = bi[i+1] - bi[i];
1475: for (k=0; k < nzL; k++) {
1476: pc1 = rtmp1 + row;
1477: pc2 = rtmp2 + row;
1478: pc3 = rtmp3 + row;
1479: if (*pc1 != 0.0 || *pc2 != 0.0 || *pc3 != 0.0) {
1480: pv = b->a + bdiag[row];
1481: mul1 = *pc1*(*pv); mul2 = *pc2*(*pv); mul3 = *pc3*(*pv);
1482: *pc1 = mul1; *pc2 = mul2; *pc3 = mul3;
1484: pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
1485: pv = b->a + bdiag[row+1]+1;
1486: nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */
1487: for (j=0; j<nz; j++) {
1488: col = pj[j];
1489: rtmp1[col] -= mul1 * pv[j];
1490: rtmp2[col] -= mul2 * pv[j];
1491: rtmp3[col] -= mul3 * pv[j];
1492: }
1493: PetscLogFlops(3+6*nz);
1494: }
1495: row = *bjtmp++;
1496: }
1498: /* finished row i; check zero pivot, then stick row i into b->a */
1499: rs = 0.0;
1500: /* L part */
1501: pc1 = b->a + bi[i];
1502: pj = b->j + bi[i];
1503: nz = bi[i+1] - bi[i];
1504: for (j=0; j<nz; j++) {
1505: col = pj[j];
1506: pc1[j] = rtmp1[col]; rs += PetscAbsScalar(pc1[j]);
1507: }
1508: /* U part */
1509: pc1 = b->a + bdiag[i+1]+1;
1510: pj = b->j + bdiag[i+1]+1;
1511: nz = bdiag[i] - bdiag[i+1] - 1; /* exclude diagonal */
1512: for (j=0; j<nz; j++) {
1513: col = pj[j];
1514: pc1[j] = rtmp1[col]; rs += PetscAbsScalar(pc1[j]);
1515: }
1517: sctx.rs = rs;
1518: sctx.pv = rtmp1[i];
1519: MatPivotCheck(B,A,info,&sctx,i);
1520: if (sctx.newshift) break;
1521: pc1 = b->a + bdiag[i]; /* Mark diag[i] */
1522: *pc1 = 1.0/sctx.pv;
1524: /* Now take care of 1st column of diagonal 3x3 block. */
1525: pc2 = rtmp2 + i;
1526: pc3 = rtmp3 + i;
1527: if (*pc2 != 0.0 || *pc3 != 0.0) {
1528: mul2 = (*pc2)*(*pc1); *pc2 = mul2;
1529: mul3 = (*pc3)*(*pc1); *pc3 = mul3;
1530: pj = b->j + bdiag[i+1]+1; /* beginning of U(i,:) */
1531: nz = bdiag[i]-bdiag[i+1]-1; /* num of entries in U(i,:) excluding diag */
1532: for (j=0; j<nz; j++) {
1533: col = pj[j];
1534: rtmp2[col] -= mul2 * rtmp1[col];
1535: rtmp3[col] -= mul3 * rtmp1[col];
1536: }
1537: PetscLogFlops(2+4*nz);
1538: }
1540: /* finished row i+1; check zero pivot, then stick row i+1 into b->a */
1541: rs = 0.0;
1542: /* L part */
1543: pc2 = b->a + bi[i+1];
1544: pj = b->j + bi[i+1];
1545: nz = bi[i+2] - bi[i+1];
1546: for (j=0; j<nz; j++) {
1547: col = pj[j];
1548: pc2[j] = rtmp2[col]; rs += PetscAbsScalar(pc2[j]);
1549: }
1550: /* U part */
1551: pc2 = b->a + bdiag[i+2]+1;
1552: pj = b->j + bdiag[i+2]+1;
1553: nz = bdiag[i+1] - bdiag[i+2] - 1; /* exclude diagonal */
1554: for (j=0; j<nz; j++) {
1555: col = pj[j];
1556: pc2[j] = rtmp2[col]; rs += PetscAbsScalar(pc2[j]);
1557: }
1559: sctx.rs = rs;
1560: sctx.pv = rtmp2[i+1];
1561: MatPivotCheck(B,A,info,&sctx,i+1);
1562: if (sctx.newshift) break;
1563: pc2 = b->a + bdiag[i+1];
1564: *pc2 = 1.0/sctx.pv; /* Mark diag[i+1] */
1566: /* Now take care of 2nd column of diagonal 3x3 block. */
1567: pc3 = rtmp3 + i+1;
1568: if (*pc3 != 0.0) {
1569: mul3 = (*pc3)*(*pc2); *pc3 = mul3;
1570: pj = b->j + bdiag[i+2]+1; /* beginning of U(i+1,:) */
1571: nz = bdiag[i+1]-bdiag[i+2]-1; /* num of entries in U(i+1,:) excluding diag */
1572: for (j=0; j<nz; j++) {
1573: col = pj[j];
1574: rtmp3[col] -= mul3 * rtmp2[col];
1575: }
1576: PetscLogFlops(1+2*nz);
1577: }
1579: /* finished i+2; check zero pivot, then stick row i+2 into b->a */
1580: rs = 0.0;
1581: /* L part */
1582: pc3 = b->a + bi[i+2];
1583: pj = b->j + bi[i+2];
1584: nz = bi[i+3] - bi[i+2];
1585: for (j=0; j<nz; j++) {
1586: col = pj[j];
1587: pc3[j] = rtmp3[col]; rs += PetscAbsScalar(pc3[j]);
1588: }
1589: /* U part */
1590: pc3 = b->a + bdiag[i+3]+1;
1591: pj = b->j + bdiag[i+3]+1;
1592: nz = bdiag[i+2] - bdiag[i+3] - 1; /* exclude diagonal */
1593: for (j=0; j<nz; j++) {
1594: col = pj[j];
1595: pc3[j] = rtmp3[col]; rs += PetscAbsScalar(pc3[j]);
1596: }
1598: sctx.rs = rs;
1599: sctx.pv = rtmp3[i+2];
1600: MatPivotCheck(B,A,info,&sctx,i+2);
1601: if (sctx.newshift) break;
1602: pc3 = b->a + bdiag[i+2];
1603: *pc3 = 1.0/sctx.pv; /* Mark diag[i+2] */
1604: break;
1605: case 4:
1606: /*----------*/
1607: /* zero rtmp */
1608: /* L part */
1609: nz = bi[i+1] - bi[i];
1610: bjtmp = bj + bi[i];
1611: for (j=0; j<nz; j++) {
1612: col = bjtmp[j];
1613: rtmp1[col] = 0.0; rtmp2[col] = 0.0; rtmp3[col] = 0.0;rtmp4[col] = 0.0;
1614: }
1616: /* U part */
1617: nz = bdiag[i]-bdiag[i+1];
1618: bjtmp = bj + bdiag[i+1]+1;
1619: for (j=0; j<nz; j++) {
1620: col = bjtmp[j];
1621: rtmp1[col] = 0.0; rtmp2[col] = 0.0; rtmp3[col] = 0.0; rtmp4[col] = 0.0;
1622: }
1624: /* load in initial (unfactored row) */
1625: nz = ai[r[i]+1] - ai[r[i]];
1626: ajtmp = aj + ai[r[i]];
1627: v1 = aa + ai[r[i]]; v2 = aa + ai[r[i]+1]; v3 = aa + ai[r[i]+2]; v4 = aa + ai[r[i]+3];
1628: for (j=0; j<nz; j++) {
1629: col = ics[ajtmp[j]];
1630: rtmp1[col] = v1[j]; rtmp2[col] = v2[j]; rtmp3[col] = v3[j]; rtmp4[col] = v4[j];
1631: }
1632: /* ZeropivotApply(): shift the diagonal of the matrix */
1633: rtmp1[i] += sctx.shift_amount; rtmp2[i+1] += sctx.shift_amount; rtmp3[i+2] += sctx.shift_amount; rtmp4[i+3] += sctx.shift_amount;
1635: /* elimination */
1636: bjtmp = bj + bi[i];
1637: row = *bjtmp++; /* pivot row */
1638: nzL = bi[i+1] - bi[i];
1639: for (k=0; k < nzL; k++) {
1640: pc1 = rtmp1 + row;
1641: pc2 = rtmp2 + row;
1642: pc3 = rtmp3 + row;
1643: pc4 = rtmp4 + row;
1644: if (*pc1 != 0.0 || *pc2 != 0.0 || *pc3 != 0.0 || *pc4 != 0.0) {
1645: pv = b->a + bdiag[row];
1646: mul1 = *pc1*(*pv); mul2 = *pc2*(*pv); mul3 = *pc3*(*pv); mul4 = *pc4*(*pv);
1647: *pc1 = mul1; *pc2 = mul2; *pc3 = mul3; *pc4 = mul4;
1649: pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */
1650: pv = b->a + bdiag[row+1]+1;
1651: nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */
1652: for (j=0; j<nz; j++) {
1653: col = pj[j];
1654: rtmp1[col] -= mul1 * pv[j];
1655: rtmp2[col] -= mul2 * pv[j];
1656: rtmp3[col] -= mul3 * pv[j];
1657: rtmp4[col] -= mul4 * pv[j];
1658: }
1659: PetscLogFlops(4+8*nz);
1660: }
1661: row = *bjtmp++;
1662: }
1664: /* finished row i; check zero pivot, then stick row i into b->a */
1665: rs = 0.0;
1666: /* L part */
1667: pc1 = b->a + bi[i];
1668: pj = b->j + bi[i];
1669: nz = bi[i+1] - bi[i];
1670: for (j=0; j<nz; j++) {
1671: col = pj[j];
1672: pc1[j] = rtmp1[col]; rs += PetscAbsScalar(pc1[j]);
1673: }
1674: /* U part */
1675: pc1 = b->a + bdiag[i+1]+1;
1676: pj = b->j + bdiag[i+1]+1;
1677: nz = bdiag[i] - bdiag[i+1] - 1; /* exclude diagonal */
1678: for (j=0; j<nz; j++) {
1679: col = pj[j];
1680: pc1[j] = rtmp1[col]; rs += PetscAbsScalar(pc1[j]);
1681: }
1683: sctx.rs = rs;
1684: sctx.pv = rtmp1[i];
1685: MatPivotCheck(B,A,info,&sctx,i);
1686: if (sctx.newshift) break;
1687: pc1 = b->a + bdiag[i]; /* Mark diag[i] */
1688: *pc1 = 1.0/sctx.pv;
1690: /* Now take care of 1st column of diagonal 4x4 block. */
1691: pc2 = rtmp2 + i;
1692: pc3 = rtmp3 + i;
1693: pc4 = rtmp4 + i;
1694: if (*pc2 != 0.0 || *pc3 != 0.0 || *pc4 != 0.0) {
1695: mul2 = (*pc2)*(*pc1); *pc2 = mul2;
1696: mul3 = (*pc3)*(*pc1); *pc3 = mul3;
1697: mul4 = (*pc4)*(*pc1); *pc4 = mul4;
1698: pj = b->j + bdiag[i+1]+1; /* beginning of U(i,:) */
1699: nz = bdiag[i]-bdiag[i+1]-1; /* num of entries in U(i,:) excluding diag */
1700: for (j=0; j<nz; j++) {
1701: col = pj[j];
1702: rtmp2[col] -= mul2 * rtmp1[col];
1703: rtmp3[col] -= mul3 * rtmp1[col];
1704: rtmp4[col] -= mul4 * rtmp1[col];
1705: }
1706: PetscLogFlops(3+6*nz);
1707: }
1709: /* finished row i+1; check zero pivot, then stick row i+1 into b->a */
1710: rs = 0.0;
1711: /* L part */
1712: pc2 = b->a + bi[i+1];
1713: pj = b->j + bi[i+1];
1714: nz = bi[i+2] - bi[i+1];
1715: for (j=0; j<nz; j++) {
1716: col = pj[j];
1717: pc2[j] = rtmp2[col]; rs += PetscAbsScalar(pc2[j]);
1718: }
1719: /* U part */
1720: pc2 = b->a + bdiag[i+2]+1;
1721: pj = b->j + bdiag[i+2]+1;
1722: nz = bdiag[i+1] - bdiag[i+2] - 1; /* exclude diagonal */
1723: for (j=0; j<nz; j++) {
1724: col = pj[j];
1725: pc2[j] = rtmp2[col]; rs += PetscAbsScalar(pc2[j]);
1726: }
1728: sctx.rs = rs;
1729: sctx.pv = rtmp2[i+1];
1730: MatPivotCheck(B,A,info,&sctx,i+1);
1731: if (sctx.newshift) break;
1732: pc2 = b->a + bdiag[i+1];
1733: *pc2 = 1.0/sctx.pv; /* Mark diag[i+1] */
1735: /* Now take care of 2nd column of diagonal 4x4 block. */
1736: pc3 = rtmp3 + i+1;
1737: pc4 = rtmp4 + i+1;
1738: if (*pc3 != 0.0 || *pc4 != 0.0) {
1739: mul3 = (*pc3)*(*pc2); *pc3 = mul3;
1740: mul4 = (*pc4)*(*pc2); *pc4 = mul4;
1741: pj = b->j + bdiag[i+2]+1; /* beginning of U(i+1,:) */
1742: nz = bdiag[i+1]-bdiag[i+2]-1; /* num of entries in U(i+1,:) excluding diag */
1743: for (j=0; j<nz; j++) {
1744: col = pj[j];
1745: rtmp3[col] -= mul3 * rtmp2[col];
1746: rtmp4[col] -= mul4 * rtmp2[col];
1747: }
1748: PetscLogFlops(4*nz);
1749: }
1751: /* finished i+2; check zero pivot, then stick row i+2 into b->a */
1752: rs = 0.0;
1753: /* L part */
1754: pc3 = b->a + bi[i+2];
1755: pj = b->j + bi[i+2];
1756: nz = bi[i+3] - bi[i+2];
1757: for (j=0; j<nz; j++) {
1758: col = pj[j];
1759: pc3[j] = rtmp3[col]; rs += PetscAbsScalar(pc3[j]);
1760: }
1761: /* U part */
1762: pc3 = b->a + bdiag[i+3]+1;
1763: pj = b->j + bdiag[i+3]+1;
1764: nz = bdiag[i+2] - bdiag[i+3] - 1; /* exclude diagonal */
1765: for (j=0; j<nz; j++) {
1766: col = pj[j];
1767: pc3[j] = rtmp3[col]; rs += PetscAbsScalar(pc3[j]);
1768: }
1770: sctx.rs = rs;
1771: sctx.pv = rtmp3[i+2];
1772: MatPivotCheck(B,A,info,&sctx,i+2);
1773: if (sctx.newshift) break;
1774: pc3 = b->a + bdiag[i+2];
1775: *pc3 = 1.0/sctx.pv; /* Mark diag[i+2] */
1777: /* Now take care of 3rd column of diagonal 4x4 block. */
1778: pc4 = rtmp4 + i+2;
1779: if (*pc4 != 0.0) {
1780: mul4 = (*pc4)*(*pc3); *pc4 = mul4;
1781: pj = b->j + bdiag[i+3]+1; /* beginning of U(i+2,:) */
1782: nz = bdiag[i+2]-bdiag[i+3]-1; /* num of entries in U(i+2,:) excluding diag */
1783: for (j=0; j<nz; j++) {
1784: col = pj[j];
1785: rtmp4[col] -= mul4 * rtmp3[col];
1786: }
1787: PetscLogFlops(1+2*nz);
1788: }
1790: /* finished i+3; check zero pivot, then stick row i+3 into b->a */
1791: rs = 0.0;
1792: /* L part */
1793: pc4 = b->a + bi[i+3];
1794: pj = b->j + bi[i+3];
1795: nz = bi[i+4] - bi[i+3];
1796: for (j=0; j<nz; j++) {
1797: col = pj[j];
1798: pc4[j] = rtmp4[col]; rs += PetscAbsScalar(pc4[j]);
1799: }
1800: /* U part */
1801: pc4 = b->a + bdiag[i+4]+1;
1802: pj = b->j + bdiag[i+4]+1;
1803: nz = bdiag[i+3] - bdiag[i+4] - 1; /* exclude diagonal */
1804: for (j=0; j<nz; j++) {
1805: col = pj[j];
1806: pc4[j] = rtmp4[col]; rs += PetscAbsScalar(pc4[j]);
1807: }
1809: sctx.rs = rs;
1810: sctx.pv = rtmp4[i+3];
1811: MatPivotCheck(B,A,info,&sctx,i+3);
1812: if (sctx.newshift) break;
1813: pc4 = b->a + bdiag[i+3];
1814: *pc4 = 1.0/sctx.pv; /* Mark diag[i+3] */
1815: break;
1817: default:
1818: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Node size not yet supported \n");
1819: }
1820: if (sctx.newshift) break; /* break for (inod=0,i=0; inod<node_max; inod++) */
1821: i += nodesz; /* Update the row */
1822: }
1824: /* MatPivotRefine() */
1825: if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) {
1826: /*
1827: * if no shift in this attempt & shifting & started shifting & can refine,
1828: * then try lower shift
1829: */
1830: sctx.shift_hi = sctx.shift_fraction;
1831: sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.;
1832: sctx.shift_amount = sctx.shift_fraction * sctx.shift_top;
1833: sctx.newshift = PETSC_TRUE;
1834: sctx.nshift++;
1835: }
1836: } while (sctx.newshift);
1838: PetscFree4(rtmp1,rtmp2,rtmp3,rtmp4);
1839: PetscFree(tmp_vec2);
1840: ISRestoreIndices(isicol,&ic);
1841: ISRestoreIndices(isrow,&r);
1843: if (b->inode.size) {
1844: C->ops->solve = MatSolve_SeqAIJ_Inode;
1845: } else {
1846: C->ops->solve = MatSolve_SeqAIJ;
1847: }
1848: C->ops->solveadd = MatSolveAdd_SeqAIJ;
1849: C->ops->solvetranspose = MatSolveTranspose_SeqAIJ;
1850: C->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ;
1851: C->ops->matsolve = MatMatSolve_SeqAIJ;
1852: C->assembled = PETSC_TRUE;
1853: C->preallocated = PETSC_TRUE;
1855: PetscLogFlops(C->cmap->n);
1857: /* MatShiftView(A,info,&sctx) */
1858: if (sctx.nshift) {
1859: if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) {
1860: PetscInfo4(A,"number of shift_pd tries %D, shift_amount %g, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,(double)sctx.shift_amount,(double)sctx.shift_fraction,(double)sctx.shift_top);
1861: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
1862: PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
1863: } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) {
1864: PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);
1865: }
1866: }
1867: return(0);
1868: }
1870: PetscErrorCode MatLUFactorNumeric_SeqAIJ_Inode_inplace(Mat B,Mat A,const MatFactorInfo *info)
1871: {
1872: Mat C = B;
1873: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)C->data;
1874: IS iscol = b->col,isrow = b->row,isicol = b->icol;
1875: PetscErrorCode ierr;
1876: const PetscInt *r,*ic,*c,*ics;
1877: PetscInt n = A->rmap->n,*bi = b->i;
1878: PetscInt *bj = b->j,*nbj=b->j +1,*ajtmp,*bjtmp,nz,nz_tmp,row,prow;
1879: PetscInt i,j,idx,*bd = b->diag,node_max,nodesz;
1880: PetscInt *ai = a->i,*aj = a->j;
1881: PetscInt *ns,*tmp_vec1,*tmp_vec2,*nsmap,*pj;
1882: PetscScalar mul1,mul2,mul3,tmp;
1883: MatScalar *pc1,*pc2,*pc3,*ba = b->a,*pv,*rtmp11,*rtmp22,*rtmp33;
1884: const MatScalar *v1,*v2,*v3,*aa = a->a,*rtmp1;
1885: PetscReal rs=0.0;
1886: FactorShiftCtx sctx;
1889: sctx.shift_top = 0;
1890: sctx.nshift_max = 0;
1891: sctx.shift_lo = 0;
1892: sctx.shift_hi = 0;
1893: sctx.shift_fraction = 0;
1895: /* if both shift schemes are chosen by user, only use info->shiftpd */
1896: if (info->shifttype==(PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */
1897: sctx.shift_top = 0;
1898: for (i=0; i<n; i++) {
1899: /* calculate rs = sum(|aij|)-RealPart(aii), amt of shift needed for this row */
1900: rs = 0.0;
1901: ajtmp = aj + ai[i];
1902: rtmp1 = aa + ai[i];
1903: nz = ai[i+1] - ai[i];
1904: for (j=0; j<nz; j++) {
1905: if (*ajtmp != i) {
1906: rs += PetscAbsScalar(*rtmp1++);
1907: } else {
1908: rs -= PetscRealPart(*rtmp1++);
1909: }
1910: ajtmp++;
1911: }
1912: if (rs>sctx.shift_top) sctx.shift_top = rs;
1913: }
1914: if (sctx.shift_top == 0.0) sctx.shift_top += 1.e-12;
1915: sctx.shift_top *= 1.1;
1916: sctx.nshift_max = 5;
1917: sctx.shift_lo = 0.;
1918: sctx.shift_hi = 1.;
1919: }
1920: sctx.shift_amount = 0;
1921: sctx.nshift = 0;
1923: ISGetIndices(isrow,&r);
1924: ISGetIndices(iscol,&c);
1925: ISGetIndices(isicol,&ic);
1926: PetscCalloc3(n,&rtmp11,n,&rtmp22,n,&rtmp33);
1927: ics = ic;
1929: node_max = a->inode.node_count;
1930: ns = a->inode.size;
1931: if (!ns) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Matrix without inode information");
1933: /* If max inode size > 3, split it into two inodes.*/
1934: /* also map the inode sizes according to the ordering */
1935: PetscMalloc1(n+1,&tmp_vec1);
1936: for (i=0,j=0; i<node_max; ++i,++j) {
1937: if (ns[i]>3) {
1938: tmp_vec1[j] = ns[i]/2; /* Assuming ns[i] < =5 */
1939: ++j;
1940: tmp_vec1[j] = ns[i] - tmp_vec1[j-1];
1941: } else {
1942: tmp_vec1[j] = ns[i];
1943: }
1944: }
1945: /* Use the correct node_max */
1946: node_max = j;
1948: /* Now reorder the inode info based on mat re-ordering info */
1949: /* First create a row -> inode_size_array_index map */
1950: PetscMalloc1(n+1,&nsmap);
1951: PetscMalloc1(node_max+1,&tmp_vec2);
1952: for (i=0,row=0; i<node_max; i++) {
1953: nodesz = tmp_vec1[i];
1954: for (j=0; j<nodesz; j++,row++) {
1955: nsmap[row] = i;
1956: }
1957: }
1958: /* Using nsmap, create a reordered ns structure */
1959: for (i=0,j=0; i< node_max; i++) {
1960: nodesz = tmp_vec1[nsmap[r[j]]]; /* here the reordered row_no is in r[] */
1961: tmp_vec2[i] = nodesz;
1962: j += nodesz;
1963: }
1964: PetscFree(nsmap);
1965: PetscFree(tmp_vec1);
1966: /* Now use the correct ns */
1967: ns = tmp_vec2;
1969: do {
1970: sctx.newshift = PETSC_FALSE;
1971: /* Now loop over each block-row, and do the factorization */
1972: for (i=0,row=0; i<node_max; i++) {
1973: nodesz = ns[i];
1974: nz = bi[row+1] - bi[row];
1975: bjtmp = bj + bi[row];
1977: switch (nodesz) {
1978: case 1:
1979: for (j=0; j<nz; j++) {
1980: idx = bjtmp[j];
1981: rtmp11[idx] = 0.0;
1982: }
1984: /* load in initial (unfactored row) */
1985: idx = r[row];
1986: nz_tmp = ai[idx+1] - ai[idx];
1987: ajtmp = aj + ai[idx];
1988: v1 = aa + ai[idx];
1990: for (j=0; j<nz_tmp; j++) {
1991: idx = ics[ajtmp[j]];
1992: rtmp11[idx] = v1[j];
1993: }
1994: rtmp11[ics[r[row]]] += sctx.shift_amount;
1996: prow = *bjtmp++;
1997: while (prow < row) {
1998: pc1 = rtmp11 + prow;
1999: if (*pc1 != 0.0) {
2000: pv = ba + bd[prow];
2001: pj = nbj + bd[prow];
2002: mul1 = *pc1 * *pv++;
2003: *pc1 = mul1;
2004: nz_tmp = bi[prow+1] - bd[prow] - 1;
2005: PetscLogFlops(1+2*nz_tmp);
2006: for (j=0; j<nz_tmp; j++) {
2007: tmp = pv[j];
2008: idx = pj[j];
2009: rtmp11[idx] -= mul1 * tmp;
2010: }
2011: }
2012: prow = *bjtmp++;
2013: }
2014: pj = bj + bi[row];
2015: pc1 = ba + bi[row];
2017: sctx.pv = rtmp11[row];
2018: rtmp11[row] = 1.0/rtmp11[row]; /* invert diag */
2019: rs = 0.0;
2020: for (j=0; j<nz; j++) {
2021: idx = pj[j];
2022: pc1[j] = rtmp11[idx]; /* rtmp11 -> ba */
2023: if (idx != row) rs += PetscAbsScalar(pc1[j]);
2024: }
2025: sctx.rs = rs;
2026: MatPivotCheck(B,A,info,&sctx,row);
2027: if (sctx.newshift) goto endofwhile;
2028: break;
2030: case 2:
2031: for (j=0; j<nz; j++) {
2032: idx = bjtmp[j];
2033: rtmp11[idx] = 0.0;
2034: rtmp22[idx] = 0.0;
2035: }
2037: /* load in initial (unfactored row) */
2038: idx = r[row];
2039: nz_tmp = ai[idx+1] - ai[idx];
2040: ajtmp = aj + ai[idx];
2041: v1 = aa + ai[idx];
2042: v2 = aa + ai[idx+1];
2043: for (j=0; j<nz_tmp; j++) {
2044: idx = ics[ajtmp[j]];
2045: rtmp11[idx] = v1[j];
2046: rtmp22[idx] = v2[j];
2047: }
2048: rtmp11[ics[r[row]]] += sctx.shift_amount;
2049: rtmp22[ics[r[row+1]]] += sctx.shift_amount;
2051: prow = *bjtmp++;
2052: while (prow < row) {
2053: pc1 = rtmp11 + prow;
2054: pc2 = rtmp22 + prow;
2055: if (*pc1 != 0.0 || *pc2 != 0.0) {
2056: pv = ba + bd[prow];
2057: pj = nbj + bd[prow];
2058: mul1 = *pc1 * *pv;
2059: mul2 = *pc2 * *pv;
2060: ++pv;
2061: *pc1 = mul1;
2062: *pc2 = mul2;
2064: nz_tmp = bi[prow+1] - bd[prow] - 1;
2065: for (j=0; j<nz_tmp; j++) {
2066: tmp = pv[j];
2067: idx = pj[j];
2068: rtmp11[idx] -= mul1 * tmp;
2069: rtmp22[idx] -= mul2 * tmp;
2070: }
2071: PetscLogFlops(2+4*nz_tmp);
2072: }
2073: prow = *bjtmp++;
2074: }
2076: /* Now take care of diagonal 2x2 block. Note: prow = row here */
2077: pc1 = rtmp11 + prow;
2078: pc2 = rtmp22 + prow;
2080: sctx.pv = *pc1;
2081: pj = bj + bi[prow];
2082: rs = 0.0;
2083: for (j=0; j<nz; j++) {
2084: idx = pj[j];
2085: if (idx != prow) rs += PetscAbsScalar(rtmp11[idx]);
2086: }
2087: sctx.rs = rs;
2088: MatPivotCheck(B,A,info,&sctx,row);
2089: if (sctx.newshift) goto endofwhile;
2091: if (*pc2 != 0.0) {
2092: pj = nbj + bd[prow];
2093: mul2 = (*pc2)/(*pc1); /* since diag is not yet inverted.*/
2094: *pc2 = mul2;
2095: nz_tmp = bi[prow+1] - bd[prow] - 1;
2096: for (j=0; j<nz_tmp; j++) {
2097: idx = pj[j];
2098: tmp = rtmp11[idx];
2099: rtmp22[idx] -= mul2 * tmp;
2100: }
2101: PetscLogFlops(1+2*nz_tmp);
2102: }
2104: pj = bj + bi[row];
2105: pc1 = ba + bi[row];
2106: pc2 = ba + bi[row+1];
2108: sctx.pv = rtmp22[row+1];
2109: rs = 0.0;
2110: rtmp11[row] = 1.0/rtmp11[row];
2111: rtmp22[row+1] = 1.0/rtmp22[row+1];
2112: /* copy row entries from dense representation to sparse */
2113: for (j=0; j<nz; j++) {
2114: idx = pj[j];
2115: pc1[j] = rtmp11[idx];
2116: pc2[j] = rtmp22[idx];
2117: if (idx != row+1) rs += PetscAbsScalar(pc2[j]);
2118: }
2119: sctx.rs = rs;
2120: MatPivotCheck(B,A,info,&sctx,row+1);
2121: if (sctx.newshift) goto endofwhile;
2122: break;
2124: case 3:
2125: for (j=0; j<nz; j++) {
2126: idx = bjtmp[j];
2127: rtmp11[idx] = 0.0;
2128: rtmp22[idx] = 0.0;
2129: rtmp33[idx] = 0.0;
2130: }
2131: /* copy the nonzeros for the 3 rows from sparse representation to dense in rtmp*[] */
2132: idx = r[row];
2133: nz_tmp = ai[idx+1] - ai[idx];
2134: ajtmp = aj + ai[idx];
2135: v1 = aa + ai[idx];
2136: v2 = aa + ai[idx+1];
2137: v3 = aa + ai[idx+2];
2138: for (j=0; j<nz_tmp; j++) {
2139: idx = ics[ajtmp[j]];
2140: rtmp11[idx] = v1[j];
2141: rtmp22[idx] = v2[j];
2142: rtmp33[idx] = v3[j];
2143: }
2144: rtmp11[ics[r[row]]] += sctx.shift_amount;
2145: rtmp22[ics[r[row+1]]] += sctx.shift_amount;
2146: rtmp33[ics[r[row+2]]] += sctx.shift_amount;
2148: /* loop over all pivot row blocks above this row block */
2149: prow = *bjtmp++;
2150: while (prow < row) {
2151: pc1 = rtmp11 + prow;
2152: pc2 = rtmp22 + prow;
2153: pc3 = rtmp33 + prow;
2154: if (*pc1 != 0.0 || *pc2 != 0.0 || *pc3 !=0.0) {
2155: pv = ba + bd[prow];
2156: pj = nbj + bd[prow];
2157: mul1 = *pc1 * *pv;
2158: mul2 = *pc2 * *pv;
2159: mul3 = *pc3 * *pv;
2160: ++pv;
2161: *pc1 = mul1;
2162: *pc2 = mul2;
2163: *pc3 = mul3;
2165: nz_tmp = bi[prow+1] - bd[prow] - 1;
2166: /* update this row based on pivot row */
2167: for (j=0; j<nz_tmp; j++) {
2168: tmp = pv[j];
2169: idx = pj[j];
2170: rtmp11[idx] -= mul1 * tmp;
2171: rtmp22[idx] -= mul2 * tmp;
2172: rtmp33[idx] -= mul3 * tmp;
2173: }
2174: PetscLogFlops(3+6*nz_tmp);
2175: }
2176: prow = *bjtmp++;
2177: }
2179: /* Now take care of diagonal 3x3 block in this set of rows */
2180: /* note: prow = row here */
2181: pc1 = rtmp11 + prow;
2182: pc2 = rtmp22 + prow;
2183: pc3 = rtmp33 + prow;
2185: sctx.pv = *pc1;
2186: pj = bj + bi[prow];
2187: rs = 0.0;
2188: for (j=0; j<nz; j++) {
2189: idx = pj[j];
2190: if (idx != row) rs += PetscAbsScalar(rtmp11[idx]);
2191: }
2192: sctx.rs = rs;
2193: MatPivotCheck(B,A,info,&sctx,row);
2194: if (sctx.newshift) goto endofwhile;
2196: if (*pc2 != 0.0 || *pc3 != 0.0) {
2197: mul2 = (*pc2)/(*pc1);
2198: mul3 = (*pc3)/(*pc1);
2199: *pc2 = mul2;
2200: *pc3 = mul3;
2201: nz_tmp = bi[prow+1] - bd[prow] - 1;
2202: pj = nbj + bd[prow];
2203: for (j=0; j<nz_tmp; j++) {
2204: idx = pj[j];
2205: tmp = rtmp11[idx];
2206: rtmp22[idx] -= mul2 * tmp;
2207: rtmp33[idx] -= mul3 * tmp;
2208: }
2209: PetscLogFlops(2+4*nz_tmp);
2210: }
2211: ++prow;
2213: pc2 = rtmp22 + prow;
2214: pc3 = rtmp33 + prow;
2215: sctx.pv = *pc2;
2216: pj = bj + bi[prow];
2217: rs = 0.0;
2218: for (j=0; j<nz; j++) {
2219: idx = pj[j];
2220: if (idx != prow) rs += PetscAbsScalar(rtmp22[idx]);
2221: }
2222: sctx.rs = rs;
2223: MatPivotCheck(B,A,info,&sctx,row+1);
2224: if (sctx.newshift) goto endofwhile;
2226: if (*pc3 != 0.0) {
2227: mul3 = (*pc3)/(*pc2);
2228: *pc3 = mul3;
2229: pj = nbj + bd[prow];
2230: nz_tmp = bi[prow+1] - bd[prow] - 1;
2231: for (j=0; j<nz_tmp; j++) {
2232: idx = pj[j];
2233: tmp = rtmp22[idx];
2234: rtmp33[idx] -= mul3 * tmp;
2235: }
2236: PetscLogFlops(1+2*nz_tmp);
2237: }
2239: pj = bj + bi[row];
2240: pc1 = ba + bi[row];
2241: pc2 = ba + bi[row+1];
2242: pc3 = ba + bi[row+2];
2244: sctx.pv = rtmp33[row+2];
2245: rs = 0.0;
2246: rtmp11[row] = 1.0/rtmp11[row];
2247: rtmp22[row+1] = 1.0/rtmp22[row+1];
2248: rtmp33[row+2] = 1.0/rtmp33[row+2];
2249: /* copy row entries from dense representation to sparse */
2250: for (j=0; j<nz; j++) {
2251: idx = pj[j];
2252: pc1[j] = rtmp11[idx];
2253: pc2[j] = rtmp22[idx];
2254: pc3[j] = rtmp33[idx];
2255: if (idx != row+2) rs += PetscAbsScalar(pc3[j]);
2256: }
2258: sctx.rs = rs;
2259: MatPivotCheck(B,A,info,&sctx,row+2);
2260: if (sctx.newshift) goto endofwhile;
2261: break;
2263: default:
2264: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Node size not yet supported \n");
2265: }
2266: row += nodesz; /* Update the row */
2267: }
2268: endofwhile:;
2269: } while (sctx.newshift);
2270: PetscFree3(rtmp11,rtmp22,rtmp33);
2271: PetscFree(tmp_vec2);
2272: ISRestoreIndices(isicol,&ic);
2273: ISRestoreIndices(isrow,&r);
2274: ISRestoreIndices(iscol,&c);
2276: (B)->ops->solve = MatSolve_SeqAIJ_inplace;
2277: /* do not set solve add, since MatSolve_Inode + Add is faster */
2278: C->ops->solvetranspose = MatSolveTranspose_SeqAIJ_inplace;
2279: C->ops->solvetransposeadd = MatSolveTransposeAdd_SeqAIJ_inplace;
2280: C->assembled = PETSC_TRUE;
2281: C->preallocated = PETSC_TRUE;
2282: if (sctx.nshift) {
2283: if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) {
2284: PetscInfo4(A,"number of shift_pd tries %D, shift_amount %g, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,(double)sctx.shift_amount,(double)sctx.shift_fraction,(double)sctx.shift_top);
2285: } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
2286: PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);
2287: }
2288: }
2289: PetscLogFlops(C->cmap->n);
2290: MatSeqAIJCheckInode(C);
2291: return(0);
2292: }
2295: /* ----------------------------------------------------------- */
2296: PetscErrorCode MatSolve_SeqAIJ_Inode(Mat A,Vec bb,Vec xx)
2297: {
2298: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2299: IS iscol = a->col,isrow = a->row;
2300: PetscErrorCode ierr;
2301: const PetscInt *r,*c,*rout,*cout;
2302: PetscInt i,j,n = A->rmap->n;
2303: PetscInt node_max,row,nsz,aii,i0,i1,nz;
2304: const PetscInt *ai = a->i,*a_j = a->j,*ns,*vi,*ad,*aj;
2305: PetscScalar *x,*tmp,*tmps,tmp0,tmp1;
2306: PetscScalar sum1,sum2,sum3,sum4,sum5;
2307: const MatScalar *v1,*v2,*v3,*v4,*v5,*a_a = a->a,*aa;
2308: const PetscScalar *b;
2311: if (!a->inode.size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Missing Inode Structure");
2312: node_max = a->inode.node_count;
2313: ns = a->inode.size; /* Node Size array */
2315: VecGetArrayRead(bb,&b);
2316: VecGetArray(xx,&x);
2317: tmp = a->solve_work;
2319: ISGetIndices(isrow,&rout); r = rout;
2320: ISGetIndices(iscol,&cout); c = cout;
2322: /* forward solve the lower triangular */
2323: tmps = tmp;
2324: aa = a_a;
2325: aj = a_j;
2326: ad = a->diag;
2328: for (i = 0,row = 0; i< node_max; ++i) {
2329: nsz = ns[i];
2330: aii = ai[row];
2331: v1 = aa + aii;
2332: vi = aj + aii;
2333: nz = ai[row+1]- ai[row];
2335: if (i < node_max-1) {
2336: /* Prefetch the indices for the next block */
2337: PetscPrefetchBlock(aj+ai[row+nsz],ai[row+nsz+1]-ai[row+nsz],0,PETSC_PREFETCH_HINT_NTA); /* indices */
2338: /* Prefetch the data for the next block */
2339: PetscPrefetchBlock(aa+ai[row+nsz],ai[row+nsz+ns[i+1]]-ai[row+nsz],0,PETSC_PREFETCH_HINT_NTA);
2340: }
2342: switch (nsz) { /* Each loop in 'case' is unrolled */
2343: case 1:
2344: sum1 = b[r[row]];
2345: for (j=0; j<nz-1; j+=2) {
2346: i0 = vi[j];
2347: i1 = vi[j+1];
2348: tmp0 = tmps[i0];
2349: tmp1 = tmps[i1];
2350: sum1 -= v1[j]*tmp0 + v1[j+1]*tmp1;
2351: }
2352: if (j == nz-1) {
2353: tmp0 = tmps[vi[j]];
2354: sum1 -= v1[j]*tmp0;
2355: }
2356: tmp[row++]=sum1;
2357: break;
2358: case 2:
2359: sum1 = b[r[row]];
2360: sum2 = b[r[row+1]];
2361: v2 = aa + ai[row+1];
2363: for (j=0; j<nz-1; j+=2) {
2364: i0 = vi[j];
2365: i1 = vi[j+1];
2366: tmp0 = tmps[i0];
2367: tmp1 = tmps[i1];
2368: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2369: sum2 -= v2[j] * tmp0 + v2[j+1] * tmp1;
2370: }
2371: if (j == nz-1) {
2372: tmp0 = tmps[vi[j]];
2373: sum1 -= v1[j] *tmp0;
2374: sum2 -= v2[j] *tmp0;
2375: }
2376: sum2 -= v2[nz] * sum1;
2377: tmp[row++]=sum1;
2378: tmp[row++]=sum2;
2379: break;
2380: case 3:
2381: sum1 = b[r[row]];
2382: sum2 = b[r[row+1]];
2383: sum3 = b[r[row+2]];
2384: v2 = aa + ai[row+1];
2385: v3 = aa + ai[row+2];
2387: for (j=0; j<nz-1; j+=2) {
2388: i0 = vi[j];
2389: i1 = vi[j+1];
2390: tmp0 = tmps[i0];
2391: tmp1 = tmps[i1];
2392: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2393: sum2 -= v2[j] * tmp0 + v2[j+1] * tmp1;
2394: sum3 -= v3[j] * tmp0 + v3[j+1] * tmp1;
2395: }
2396: if (j == nz-1) {
2397: tmp0 = tmps[vi[j]];
2398: sum1 -= v1[j] *tmp0;
2399: sum2 -= v2[j] *tmp0;
2400: sum3 -= v3[j] *tmp0;
2401: }
2402: sum2 -= v2[nz] * sum1;
2403: sum3 -= v3[nz] * sum1;
2404: sum3 -= v3[nz+1] * sum2;
2405: tmp[row++]=sum1;
2406: tmp[row++]=sum2;
2407: tmp[row++]=sum3;
2408: break;
2410: case 4:
2411: sum1 = b[r[row]];
2412: sum2 = b[r[row+1]];
2413: sum3 = b[r[row+2]];
2414: sum4 = b[r[row+3]];
2415: v2 = aa + ai[row+1];
2416: v3 = aa + ai[row+2];
2417: v4 = aa + ai[row+3];
2419: for (j=0; j<nz-1; j+=2) {
2420: i0 = vi[j];
2421: i1 = vi[j+1];
2422: tmp0 = tmps[i0];
2423: tmp1 = tmps[i1];
2424: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2425: sum2 -= v2[j] * tmp0 + v2[j+1] * tmp1;
2426: sum3 -= v3[j] * tmp0 + v3[j+1] * tmp1;
2427: sum4 -= v4[j] * tmp0 + v4[j+1] * tmp1;
2428: }
2429: if (j == nz-1) {
2430: tmp0 = tmps[vi[j]];
2431: sum1 -= v1[j] *tmp0;
2432: sum2 -= v2[j] *tmp0;
2433: sum3 -= v3[j] *tmp0;
2434: sum4 -= v4[j] *tmp0;
2435: }
2436: sum2 -= v2[nz] * sum1;
2437: sum3 -= v3[nz] * sum1;
2438: sum4 -= v4[nz] * sum1;
2439: sum3 -= v3[nz+1] * sum2;
2440: sum4 -= v4[nz+1] * sum2;
2441: sum4 -= v4[nz+2] * sum3;
2443: tmp[row++]=sum1;
2444: tmp[row++]=sum2;
2445: tmp[row++]=sum3;
2446: tmp[row++]=sum4;
2447: break;
2448: case 5:
2449: sum1 = b[r[row]];
2450: sum2 = b[r[row+1]];
2451: sum3 = b[r[row+2]];
2452: sum4 = b[r[row+3]];
2453: sum5 = b[r[row+4]];
2454: v2 = aa + ai[row+1];
2455: v3 = aa + ai[row+2];
2456: v4 = aa + ai[row+3];
2457: v5 = aa + ai[row+4];
2459: for (j=0; j<nz-1; j+=2) {
2460: i0 = vi[j];
2461: i1 = vi[j+1];
2462: tmp0 = tmps[i0];
2463: tmp1 = tmps[i1];
2464: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2465: sum2 -= v2[j] * tmp0 + v2[j+1] * tmp1;
2466: sum3 -= v3[j] * tmp0 + v3[j+1] * tmp1;
2467: sum4 -= v4[j] * tmp0 + v4[j+1] * tmp1;
2468: sum5 -= v5[j] * tmp0 + v5[j+1] * tmp1;
2469: }
2470: if (j == nz-1) {
2471: tmp0 = tmps[vi[j]];
2472: sum1 -= v1[j] *tmp0;
2473: sum2 -= v2[j] *tmp0;
2474: sum3 -= v3[j] *tmp0;
2475: sum4 -= v4[j] *tmp0;
2476: sum5 -= v5[j] *tmp0;
2477: }
2479: sum2 -= v2[nz] * sum1;
2480: sum3 -= v3[nz] * sum1;
2481: sum4 -= v4[nz] * sum1;
2482: sum5 -= v5[nz] * sum1;
2483: sum3 -= v3[nz+1] * sum2;
2484: sum4 -= v4[nz+1] * sum2;
2485: sum5 -= v5[nz+1] * sum2;
2486: sum4 -= v4[nz+2] * sum3;
2487: sum5 -= v5[nz+2] * sum3;
2488: sum5 -= v5[nz+3] * sum4;
2490: tmp[row++]=sum1;
2491: tmp[row++]=sum2;
2492: tmp[row++]=sum3;
2493: tmp[row++]=sum4;
2494: tmp[row++]=sum5;
2495: break;
2496: default:
2497: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Node size not yet supported \n");
2498: }
2499: }
2500: /* backward solve the upper triangular */
2501: for (i=node_max -1,row = n-1; i>=0; i--) {
2502: nsz = ns[i];
2503: aii = ad[row+1] + 1;
2504: v1 = aa + aii;
2505: vi = aj + aii;
2506: nz = ad[row]- ad[row+1] - 1;
2508: if (i > 0) {
2509: /* Prefetch the indices for the next block */
2510: PetscPrefetchBlock(aj+ad[row-nsz+1]+1,ad[row-nsz]-ad[row-nsz+1],0,PETSC_PREFETCH_HINT_NTA);
2511: /* Prefetch the data for the next block */
2512: PetscPrefetchBlock(aa+ad[row-nsz+1]+1,ad[row-nsz-ns[i-1]+1]-ad[row-nsz+1],0,PETSC_PREFETCH_HINT_NTA);
2513: }
2515: switch (nsz) { /* Each loop in 'case' is unrolled */
2516: case 1:
2517: sum1 = tmp[row];
2519: for (j=0; j<nz-1; j+=2) {
2520: i0 = vi[j];
2521: i1 = vi[j+1];
2522: tmp0 = tmps[i0];
2523: tmp1 = tmps[i1];
2524: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2525: }
2526: if (j == nz-1) {
2527: tmp0 = tmps[vi[j]];
2528: sum1 -= v1[j]*tmp0;
2529: }
2530: x[c[row]] = tmp[row] = sum1*v1[nz]; row--;
2531: break;
2532: case 2:
2533: sum1 = tmp[row];
2534: sum2 = tmp[row-1];
2535: v2 = aa + ad[row] + 1;
2536: for (j=0; j<nz-1; j+=2) {
2537: i0 = vi[j];
2538: i1 = vi[j+1];
2539: tmp0 = tmps[i0];
2540: tmp1 = tmps[i1];
2541: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2542: sum2 -= v2[j+1] * tmp0 + v2[j+2] * tmp1;
2543: }
2544: if (j == nz-1) {
2545: tmp0 = tmps[vi[j]];
2546: sum1 -= v1[j]* tmp0;
2547: sum2 -= v2[j+1]* tmp0;
2548: }
2550: tmp0 = x[c[row]] = tmp[row] = sum1*v1[nz]; row--;
2551: sum2 -= v2[0] * tmp0;
2552: x[c[row]] = tmp[row] = sum2*v2[nz+1]; row--;
2553: break;
2554: case 3:
2555: sum1 = tmp[row];
2556: sum2 = tmp[row -1];
2557: sum3 = tmp[row -2];
2558: v2 = aa + ad[row] + 1;
2559: v3 = aa + ad[row -1] + 1;
2560: for (j=0; j<nz-1; j+=2) {
2561: i0 = vi[j];
2562: i1 = vi[j+1];
2563: tmp0 = tmps[i0];
2564: tmp1 = tmps[i1];
2565: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2566: sum2 -= v2[j+1] * tmp0 + v2[j+2] * tmp1;
2567: sum3 -= v3[j+2] * tmp0 + v3[j+3] * tmp1;
2568: }
2569: if (j== nz-1) {
2570: tmp0 = tmps[vi[j]];
2571: sum1 -= v1[j] * tmp0;
2572: sum2 -= v2[j+1] * tmp0;
2573: sum3 -= v3[j+2] * tmp0;
2574: }
2575: tmp0 = x[c[row]] = tmp[row] = sum1*v1[nz]; row--;
2576: sum2 -= v2[0]* tmp0;
2577: sum3 -= v3[1] * tmp0;
2578: tmp0 = x[c[row]] = tmp[row] = sum2*v2[nz+1]; row--;
2579: sum3 -= v3[0]* tmp0;
2580: x[c[row]] = tmp[row] = sum3*v3[nz+2]; row--;
2582: break;
2583: case 4:
2584: sum1 = tmp[row];
2585: sum2 = tmp[row -1];
2586: sum3 = tmp[row -2];
2587: sum4 = tmp[row -3];
2588: v2 = aa + ad[row]+1;
2589: v3 = aa + ad[row -1]+1;
2590: v4 = aa + ad[row -2]+1;
2592: for (j=0; j<nz-1; j+=2) {
2593: i0 = vi[j];
2594: i1 = vi[j+1];
2595: tmp0 = tmps[i0];
2596: tmp1 = tmps[i1];
2597: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2598: sum2 -= v2[j+1] * tmp0 + v2[j+2] * tmp1;
2599: sum3 -= v3[j+2] * tmp0 + v3[j+3] * tmp1;
2600: sum4 -= v4[j+3] * tmp0 + v4[j+4] * tmp1;
2601: }
2602: if (j== nz-1) {
2603: tmp0 = tmps[vi[j]];
2604: sum1 -= v1[j] * tmp0;
2605: sum2 -= v2[j+1] * tmp0;
2606: sum3 -= v3[j+2] * tmp0;
2607: sum4 -= v4[j+3] * tmp0;
2608: }
2610: tmp0 = x[c[row]] = tmp[row] = sum1*v1[nz]; row--;
2611: sum2 -= v2[0] * tmp0;
2612: sum3 -= v3[1] * tmp0;
2613: sum4 -= v4[2] * tmp0;
2614: tmp0 = x[c[row]] = tmp[row] = sum2*v2[nz+1]; row--;
2615: sum3 -= v3[0] * tmp0;
2616: sum4 -= v4[1] * tmp0;
2617: tmp0 = x[c[row]] = tmp[row] = sum3*v3[nz+2]; row--;
2618: sum4 -= v4[0] * tmp0;
2619: x[c[row]] = tmp[row] = sum4*v4[nz+3]; row--;
2620: break;
2621: case 5:
2622: sum1 = tmp[row];
2623: sum2 = tmp[row -1];
2624: sum3 = tmp[row -2];
2625: sum4 = tmp[row -3];
2626: sum5 = tmp[row -4];
2627: v2 = aa + ad[row]+1;
2628: v3 = aa + ad[row -1]+1;
2629: v4 = aa + ad[row -2]+1;
2630: v5 = aa + ad[row -3]+1;
2631: for (j=0; j<nz-1; j+=2) {
2632: i0 = vi[j];
2633: i1 = vi[j+1];
2634: tmp0 = tmps[i0];
2635: tmp1 = tmps[i1];
2636: sum1 -= v1[j] * tmp0 + v1[j+1] * tmp1;
2637: sum2 -= v2[j+1] * tmp0 + v2[j+2] * tmp1;
2638: sum3 -= v3[j+2] * tmp0 + v3[j+3] * tmp1;
2639: sum4 -= v4[j+3] * tmp0 + v4[j+4] * tmp1;
2640: sum5 -= v5[j+4] * tmp0 + v5[j+5] * tmp1;
2641: }
2642: if (j==nz-1) {
2643: tmp0 = tmps[vi[j]];
2644: sum1 -= v1[j] * tmp0;
2645: sum2 -= v2[j+1] * tmp0;
2646: sum3 -= v3[j+2] * tmp0;
2647: sum4 -= v4[j+3] * tmp0;
2648: sum5 -= v5[j+4] * tmp0;
2649: }
2651: tmp0 = x[c[row]] = tmp[row] = sum1*v1[nz]; row--;
2652: sum2 -= v2[0] * tmp0;
2653: sum3 -= v3[1] * tmp0;
2654: sum4 -= v4[2] * tmp0;
2655: sum5 -= v5[3] * tmp0;
2656: tmp0 = x[c[row]] = tmp[row] = sum2*v2[nz+1]; row--;
2657: sum3 -= v3[0] * tmp0;
2658: sum4 -= v4[1] * tmp0;
2659: sum5 -= v5[2] * tmp0;
2660: tmp0 = x[c[row]] = tmp[row] = sum3*v3[nz+2]; row--;
2661: sum4 -= v4[0] * tmp0;
2662: sum5 -= v5[1] * tmp0;
2663: tmp0 = x[c[row]] = tmp[row] = sum4*v4[nz+3]; row--;
2664: sum5 -= v5[0] * tmp0;
2665: x[c[row]] = tmp[row] = sum5*v5[nz+4]; row--;
2666: break;
2667: default:
2668: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"Node size not yet supported \n");
2669: }
2670: }
2671: ISRestoreIndices(isrow,&rout);
2672: ISRestoreIndices(iscol,&cout);
2673: VecRestoreArrayRead(bb,&b);
2674: VecRestoreArray(xx,&x);
2675: PetscLogFlops(2.0*a->nz - A->cmap->n);
2676: return(0);
2677: }
2680: /*
2681: Makes a longer coloring[] array and calls the usual code with that
2682: */
2683: PetscErrorCode MatColoringPatch_SeqAIJ_Inode(Mat mat,PetscInt ncolors,PetscInt nin,ISColoringValue coloring[],ISColoring *iscoloring)
2684: {
2685: Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data;
2686: PetscErrorCode ierr;
2687: PetscInt n = mat->cmap->n,m = a->inode.node_count,j,*ns = a->inode.size,row;
2688: PetscInt *colorused,i;
2689: ISColoringValue *newcolor;
2692: PetscMalloc1(n+1,&newcolor);
2693: /* loop over inodes, marking a color for each column*/
2694: row = 0;
2695: for (i=0; i<m; i++) {
2696: for (j=0; j<ns[i]; j++) {
2697: newcolor[row++] = coloring[i] + j*ncolors;
2698: }
2699: }
2701: /* eliminate unneeded colors */
2702: PetscCalloc1(5*ncolors,&colorused);
2703: for (i=0; i<n; i++) {
2704: colorused[newcolor[i]] = 1;
2705: }
2707: for (i=1; i<5*ncolors; i++) {
2708: colorused[i] += colorused[i-1];
2709: }
2710: ncolors = colorused[5*ncolors-1];
2711: for (i=0; i<n; i++) {
2712: newcolor[i] = colorused[newcolor[i]]-1;
2713: }
2714: PetscFree(colorused);
2715: ISColoringCreate(PetscObjectComm((PetscObject)mat),ncolors,n,newcolor,PETSC_OWN_POINTER,iscoloring);
2716: PetscFree(coloring);
2717: return(0);
2718: }
2720: #include <petsc/private/kernels/blockinvert.h>
2722: PetscErrorCode MatSOR_SeqAIJ_Inode(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
2723: {
2724: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
2725: PetscScalar sum1 = 0.0,sum2 = 0.0,sum3 = 0.0,sum4 = 0.0,sum5 = 0.0,tmp0,tmp1,tmp2,tmp3;
2726: MatScalar *ibdiag,*bdiag,work[25],*t;
2727: PetscScalar *x,tmp4,tmp5,x1,x2,x3,x4,x5;
2728: const MatScalar *v = a->a,*v1 = NULL,*v2 = NULL,*v3 = NULL,*v4 = NULL,*v5 = NULL;
2729: const PetscScalar *xb, *b;
2730: PetscReal zeropivot = 100.*PETSC_MACHINE_EPSILON, shift = 0.0;
2731: PetscErrorCode ierr;
2732: PetscInt n,m = a->inode.node_count,cnt = 0,i,j,row,i1,i2;
2733: PetscInt sz,k,ipvt[5];
2734: PetscBool allowzeropivot,zeropivotdetected;
2735: const PetscInt *sizes = a->inode.size,*idx,*diag = a->diag,*ii = a->i;
2738: allowzeropivot = PetscNot(A->erroriffailure);
2739: if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for omega != 1.0; use -mat_no_inode");
2740: if (fshift != 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for fshift != 0.0; use -mat_no_inode");
2742: if (!a->inode.ibdiagvalid) {
2743: if (!a->inode.ibdiag) {
2744: /* calculate space needed for diagonal blocks */
2745: for (i=0; i<m; i++) {
2746: cnt += sizes[i]*sizes[i];
2747: }
2748: a->inode.bdiagsize = cnt;
2750: PetscMalloc3(cnt,&a->inode.ibdiag,cnt,&a->inode.bdiag,A->rmap->n,&a->inode.ssor_work);
2751: }
2753: /* copy over the diagonal blocks and invert them */
2754: ibdiag = a->inode.ibdiag;
2755: bdiag = a->inode.bdiag;
2756: cnt = 0;
2757: for (i=0, row = 0; i<m; i++) {
2758: for (j=0; j<sizes[i]; j++) {
2759: for (k=0; k<sizes[i]; k++) {
2760: bdiag[cnt+k*sizes[i]+j] = v[diag[row+j] - j + k];
2761: }
2762: }
2763: PetscMemcpy(ibdiag+cnt,bdiag+cnt,sizes[i]*sizes[i]*sizeof(MatScalar));
2765: switch (sizes[i]) {
2766: case 1:
2767: /* Create matrix data structure */
2768: if (PetscAbsScalar(ibdiag[cnt]) < zeropivot) {
2769: if (allowzeropivot) {
2770: A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
2771: A->factorerror_zeropivot_value = PetscAbsScalar(ibdiag[cnt]);
2772: A->factorerror_zeropivot_row = row;
2773: PetscInfo1(A,"Zero pivot, row %D\n",row);
2774: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot on row %D",row);
2775: }
2776: ibdiag[cnt] = 1.0/ibdiag[cnt];
2777: break;
2778: case 2:
2779: PetscKernel_A_gets_inverse_A_2(ibdiag+cnt,shift,allowzeropivot,&zeropivotdetected);
2780: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
2781: break;
2782: case 3:
2783: PetscKernel_A_gets_inverse_A_3(ibdiag+cnt,shift,allowzeropivot,&zeropivotdetected);
2784: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
2785: break;
2786: case 4:
2787: PetscKernel_A_gets_inverse_A_4(ibdiag+cnt,shift,allowzeropivot,&zeropivotdetected);
2788: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
2789: break;
2790: case 5:
2791: PetscKernel_A_gets_inverse_A_5(ibdiag+cnt,ipvt,work,shift,allowzeropivot,&zeropivotdetected);
2792: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
2793: break;
2794: default:
2795: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
2796: }
2797: cnt += sizes[i]*sizes[i];
2798: row += sizes[i];
2799: }
2800: a->inode.ibdiagvalid = PETSC_TRUE;
2801: }
2802: ibdiag = a->inode.ibdiag;
2803: bdiag = a->inode.bdiag;
2804: t = a->inode.ssor_work;
2806: VecGetArray(xx,&x);
2807: VecGetArrayRead(bb,&b);
2808: /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */
2809: if (flag & SOR_ZERO_INITIAL_GUESS) {
2810: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
2812: for (i=0, row=0; i<m; i++) {
2813: sz = diag[row] - ii[row];
2814: v1 = a->a + ii[row];
2815: idx = a->j + ii[row];
2817: /* see comments for MatMult_SeqAIJ_Inode() for how this is coded */
2818: switch (sizes[i]) {
2819: case 1:
2821: sum1 = b[row];
2822: for (n = 0; n<sz-1; n+=2) {
2823: i1 = idx[0];
2824: i2 = idx[1];
2825: idx += 2;
2826: tmp0 = x[i1];
2827: tmp1 = x[i2];
2828: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
2829: }
2831: if (n == sz-1) {
2832: tmp0 = x[*idx];
2833: sum1 -= *v1 * tmp0;
2834: }
2835: t[row] = sum1;
2836: x[row++] = sum1*(*ibdiag++);
2837: break;
2838: case 2:
2839: v2 = a->a + ii[row+1];
2840: sum1 = b[row];
2841: sum2 = b[row+1];
2842: for (n = 0; n<sz-1; n+=2) {
2843: i1 = idx[0];
2844: i2 = idx[1];
2845: idx += 2;
2846: tmp0 = x[i1];
2847: tmp1 = x[i2];
2848: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
2849: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
2850: }
2852: if (n == sz-1) {
2853: tmp0 = x[*idx];
2854: sum1 -= v1[0] * tmp0;
2855: sum2 -= v2[0] * tmp0;
2856: }
2857: t[row] = sum1;
2858: t[row+1] = sum2;
2859: x[row++] = sum1*ibdiag[0] + sum2*ibdiag[2];
2860: x[row++] = sum1*ibdiag[1] + sum2*ibdiag[3];
2861: ibdiag += 4;
2862: break;
2863: case 3:
2864: v2 = a->a + ii[row+1];
2865: v3 = a->a + ii[row+2];
2866: sum1 = b[row];
2867: sum2 = b[row+1];
2868: sum3 = b[row+2];
2869: for (n = 0; n<sz-1; n+=2) {
2870: i1 = idx[0];
2871: i2 = idx[1];
2872: idx += 2;
2873: tmp0 = x[i1];
2874: tmp1 = x[i2];
2875: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
2876: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
2877: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
2878: }
2880: if (n == sz-1) {
2881: tmp0 = x[*idx];
2882: sum1 -= v1[0] * tmp0;
2883: sum2 -= v2[0] * tmp0;
2884: sum3 -= v3[0] * tmp0;
2885: }
2886: t[row] = sum1;
2887: t[row+1] = sum2;
2888: t[row+2] = sum3;
2889: x[row++] = sum1*ibdiag[0] + sum2*ibdiag[3] + sum3*ibdiag[6];
2890: x[row++] = sum1*ibdiag[1] + sum2*ibdiag[4] + sum3*ibdiag[7];
2891: x[row++] = sum1*ibdiag[2] + sum2*ibdiag[5] + sum3*ibdiag[8];
2892: ibdiag += 9;
2893: break;
2894: case 4:
2895: v2 = a->a + ii[row+1];
2896: v3 = a->a + ii[row+2];
2897: v4 = a->a + ii[row+3];
2898: sum1 = b[row];
2899: sum2 = b[row+1];
2900: sum3 = b[row+2];
2901: sum4 = b[row+3];
2902: for (n = 0; n<sz-1; n+=2) {
2903: i1 = idx[0];
2904: i2 = idx[1];
2905: idx += 2;
2906: tmp0 = x[i1];
2907: tmp1 = x[i2];
2908: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
2909: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
2910: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
2911: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
2912: }
2914: if (n == sz-1) {
2915: tmp0 = x[*idx];
2916: sum1 -= v1[0] * tmp0;
2917: sum2 -= v2[0] * tmp0;
2918: sum3 -= v3[0] * tmp0;
2919: sum4 -= v4[0] * tmp0;
2920: }
2921: t[row] = sum1;
2922: t[row+1] = sum2;
2923: t[row+2] = sum3;
2924: t[row+3] = sum4;
2925: x[row++] = sum1*ibdiag[0] + sum2*ibdiag[4] + sum3*ibdiag[8] + sum4*ibdiag[12];
2926: x[row++] = sum1*ibdiag[1] + sum2*ibdiag[5] + sum3*ibdiag[9] + sum4*ibdiag[13];
2927: x[row++] = sum1*ibdiag[2] + sum2*ibdiag[6] + sum3*ibdiag[10] + sum4*ibdiag[14];
2928: x[row++] = sum1*ibdiag[3] + sum2*ibdiag[7] + sum3*ibdiag[11] + sum4*ibdiag[15];
2929: ibdiag += 16;
2930: break;
2931: case 5:
2932: v2 = a->a + ii[row+1];
2933: v3 = a->a + ii[row+2];
2934: v4 = a->a + ii[row+3];
2935: v5 = a->a + ii[row+4];
2936: sum1 = b[row];
2937: sum2 = b[row+1];
2938: sum3 = b[row+2];
2939: sum4 = b[row+3];
2940: sum5 = b[row+4];
2941: for (n = 0; n<sz-1; n+=2) {
2942: i1 = idx[0];
2943: i2 = idx[1];
2944: idx += 2;
2945: tmp0 = x[i1];
2946: tmp1 = x[i2];
2947: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
2948: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
2949: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
2950: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
2951: sum5 -= v5[0] * tmp0 + v5[1] * tmp1; v5 += 2;
2952: }
2954: if (n == sz-1) {
2955: tmp0 = x[*idx];
2956: sum1 -= v1[0] * tmp0;
2957: sum2 -= v2[0] * tmp0;
2958: sum3 -= v3[0] * tmp0;
2959: sum4 -= v4[0] * tmp0;
2960: sum5 -= v5[0] * tmp0;
2961: }
2962: t[row] = sum1;
2963: t[row+1] = sum2;
2964: t[row+2] = sum3;
2965: t[row+3] = sum4;
2966: t[row+4] = sum5;
2967: x[row++] = sum1*ibdiag[0] + sum2*ibdiag[5] + sum3*ibdiag[10] + sum4*ibdiag[15] + sum5*ibdiag[20];
2968: x[row++] = sum1*ibdiag[1] + sum2*ibdiag[6] + sum3*ibdiag[11] + sum4*ibdiag[16] + sum5*ibdiag[21];
2969: x[row++] = sum1*ibdiag[2] + sum2*ibdiag[7] + sum3*ibdiag[12] + sum4*ibdiag[17] + sum5*ibdiag[22];
2970: x[row++] = sum1*ibdiag[3] + sum2*ibdiag[8] + sum3*ibdiag[13] + sum4*ibdiag[18] + sum5*ibdiag[23];
2971: x[row++] = sum1*ibdiag[4] + sum2*ibdiag[9] + sum3*ibdiag[14] + sum4*ibdiag[19] + sum5*ibdiag[24];
2972: ibdiag += 25;
2973: break;
2974: default:
2975: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
2976: }
2977: }
2979: xb = t;
2980: PetscLogFlops(a->nz);
2981: } else xb = b;
2982: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
2984: ibdiag = a->inode.ibdiag+a->inode.bdiagsize;
2985: for (i=m-1, row=A->rmap->n-1; i>=0; i--) {
2986: ibdiag -= sizes[i]*sizes[i];
2987: sz = ii[row+1] - diag[row] - 1;
2988: v1 = a->a + diag[row] + 1;
2989: idx = a->j + diag[row] + 1;
2991: /* see comments for MatMult_SeqAIJ_Inode() for how this is coded */
2992: switch (sizes[i]) {
2993: case 1:
2995: sum1 = xb[row];
2996: for (n = 0; n<sz-1; n+=2) {
2997: i1 = idx[0];
2998: i2 = idx[1];
2999: idx += 2;
3000: tmp0 = x[i1];
3001: tmp1 = x[i2];
3002: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3003: }
3005: if (n == sz-1) {
3006: tmp0 = x[*idx];
3007: sum1 -= *v1*tmp0;
3008: }
3009: x[row--] = sum1*(*ibdiag);
3010: break;
3012: case 2:
3014: sum1 = xb[row];
3015: sum2 = xb[row-1];
3016: /* note that sum1 is associated with the second of the two rows */
3017: v2 = a->a + diag[row-1] + 2;
3018: for (n = 0; n<sz-1; n+=2) {
3019: i1 = idx[0];
3020: i2 = idx[1];
3021: idx += 2;
3022: tmp0 = x[i1];
3023: tmp1 = x[i2];
3024: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3025: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3026: }
3028: if (n == sz-1) {
3029: tmp0 = x[*idx];
3030: sum1 -= *v1*tmp0;
3031: sum2 -= *v2*tmp0;
3032: }
3033: x[row--] = sum2*ibdiag[1] + sum1*ibdiag[3];
3034: x[row--] = sum2*ibdiag[0] + sum1*ibdiag[2];
3035: break;
3036: case 3:
3038: sum1 = xb[row];
3039: sum2 = xb[row-1];
3040: sum3 = xb[row-2];
3041: v2 = a->a + diag[row-1] + 2;
3042: v3 = a->a + diag[row-2] + 3;
3043: for (n = 0; n<sz-1; n+=2) {
3044: i1 = idx[0];
3045: i2 = idx[1];
3046: idx += 2;
3047: tmp0 = x[i1];
3048: tmp1 = x[i2];
3049: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3050: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3051: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3052: }
3054: if (n == sz-1) {
3055: tmp0 = x[*idx];
3056: sum1 -= *v1*tmp0;
3057: sum2 -= *v2*tmp0;
3058: sum3 -= *v3*tmp0;
3059: }
3060: x[row--] = sum3*ibdiag[2] + sum2*ibdiag[5] + sum1*ibdiag[8];
3061: x[row--] = sum3*ibdiag[1] + sum2*ibdiag[4] + sum1*ibdiag[7];
3062: x[row--] = sum3*ibdiag[0] + sum2*ibdiag[3] + sum1*ibdiag[6];
3063: break;
3064: case 4:
3066: sum1 = xb[row];
3067: sum2 = xb[row-1];
3068: sum3 = xb[row-2];
3069: sum4 = xb[row-3];
3070: v2 = a->a + diag[row-1] + 2;
3071: v3 = a->a + diag[row-2] + 3;
3072: v4 = a->a + diag[row-3] + 4;
3073: for (n = 0; n<sz-1; n+=2) {
3074: i1 = idx[0];
3075: i2 = idx[1];
3076: idx += 2;
3077: tmp0 = x[i1];
3078: tmp1 = x[i2];
3079: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3080: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3081: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3082: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3083: }
3085: if (n == sz-1) {
3086: tmp0 = x[*idx];
3087: sum1 -= *v1*tmp0;
3088: sum2 -= *v2*tmp0;
3089: sum3 -= *v3*tmp0;
3090: sum4 -= *v4*tmp0;
3091: }
3092: x[row--] = sum4*ibdiag[3] + sum3*ibdiag[7] + sum2*ibdiag[11] + sum1*ibdiag[15];
3093: x[row--] = sum4*ibdiag[2] + sum3*ibdiag[6] + sum2*ibdiag[10] + sum1*ibdiag[14];
3094: x[row--] = sum4*ibdiag[1] + sum3*ibdiag[5] + sum2*ibdiag[9] + sum1*ibdiag[13];
3095: x[row--] = sum4*ibdiag[0] + sum3*ibdiag[4] + sum2*ibdiag[8] + sum1*ibdiag[12];
3096: break;
3097: case 5:
3099: sum1 = xb[row];
3100: sum2 = xb[row-1];
3101: sum3 = xb[row-2];
3102: sum4 = xb[row-3];
3103: sum5 = xb[row-4];
3104: v2 = a->a + diag[row-1] + 2;
3105: v3 = a->a + diag[row-2] + 3;
3106: v4 = a->a + diag[row-3] + 4;
3107: v5 = a->a + diag[row-4] + 5;
3108: for (n = 0; n<sz-1; n+=2) {
3109: i1 = idx[0];
3110: i2 = idx[1];
3111: idx += 2;
3112: tmp0 = x[i1];
3113: tmp1 = x[i2];
3114: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3115: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3116: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3117: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3118: sum5 -= v5[0] * tmp0 + v5[1] * tmp1; v5 += 2;
3119: }
3121: if (n == sz-1) {
3122: tmp0 = x[*idx];
3123: sum1 -= *v1*tmp0;
3124: sum2 -= *v2*tmp0;
3125: sum3 -= *v3*tmp0;
3126: sum4 -= *v4*tmp0;
3127: sum5 -= *v5*tmp0;
3128: }
3129: x[row--] = sum5*ibdiag[4] + sum4*ibdiag[9] + sum3*ibdiag[14] + sum2*ibdiag[19] + sum1*ibdiag[24];
3130: x[row--] = sum5*ibdiag[3] + sum4*ibdiag[8] + sum3*ibdiag[13] + sum2*ibdiag[18] + sum1*ibdiag[23];
3131: x[row--] = sum5*ibdiag[2] + sum4*ibdiag[7] + sum3*ibdiag[12] + sum2*ibdiag[17] + sum1*ibdiag[22];
3132: x[row--] = sum5*ibdiag[1] + sum4*ibdiag[6] + sum3*ibdiag[11] + sum2*ibdiag[16] + sum1*ibdiag[21];
3133: x[row--] = sum5*ibdiag[0] + sum4*ibdiag[5] + sum3*ibdiag[10] + sum2*ibdiag[15] + sum1*ibdiag[20];
3134: break;
3135: default:
3136: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
3137: }
3138: }
3140: PetscLogFlops(a->nz);
3141: }
3142: its--;
3143: }
3144: while (its--) {
3146: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
3147: for (i=0, row=0, ibdiag = a->inode.ibdiag;
3148: i<m;
3149: row += sizes[i], ibdiag += sizes[i]*sizes[i], i++) {
3151: sz = diag[row] - ii[row];
3152: v1 = a->a + ii[row];
3153: idx = a->j + ii[row];
3154: /* see comments for MatMult_SeqAIJ_Inode() for how this is coded */
3155: switch (sizes[i]) {
3156: case 1:
3157: sum1 = b[row];
3158: for (n = 0; n<sz-1; n+=2) {
3159: i1 = idx[0];
3160: i2 = idx[1];
3161: idx += 2;
3162: tmp0 = x[i1];
3163: tmp1 = x[i2];
3164: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3165: }
3166: if (n == sz-1) {
3167: tmp0 = x[*idx++];
3168: sum1 -= *v1 * tmp0;
3169: v1++;
3170: }
3171: t[row] = sum1;
3172: sz = ii[row+1] - diag[row] - 1;
3173: idx = a->j + diag[row] + 1;
3174: v1 += 1;
3175: for (n = 0; n<sz-1; n+=2) {
3176: i1 = idx[0];
3177: i2 = idx[1];
3178: idx += 2;
3179: tmp0 = x[i1];
3180: tmp1 = x[i2];
3181: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3182: }
3183: if (n == sz-1) {
3184: tmp0 = x[*idx++];
3185: sum1 -= *v1 * tmp0;
3186: }
3187: /* in MatSOR_SeqAIJ this line would be
3188: *
3189: * x[row] = (1-omega)*x[row]+(sum1+(*bdiag++)*x[row])*(*ibdiag++);
3190: *
3191: * but omega == 1, so this becomes
3192: *
3193: * x[row] = sum1*(*ibdiag++);
3194: *
3195: */
3196: x[row] = sum1*(*ibdiag);
3197: break;
3198: case 2:
3199: v2 = a->a + ii[row+1];
3200: sum1 = b[row];
3201: sum2 = b[row+1];
3202: for (n = 0; n<sz-1; n+=2) {
3203: i1 = idx[0];
3204: i2 = idx[1];
3205: idx += 2;
3206: tmp0 = x[i1];
3207: tmp1 = x[i2];
3208: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3209: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3210: }
3211: if (n == sz-1) {
3212: tmp0 = x[*idx++];
3213: sum1 -= v1[0] * tmp0;
3214: sum2 -= v2[0] * tmp0;
3215: v1++; v2++;
3216: }
3217: t[row] = sum1;
3218: t[row+1] = sum2;
3219: sz = ii[row+1] - diag[row] - 2;
3220: idx = a->j + diag[row] + 2;
3221: v1 += 2;
3222: v2 += 2;
3223: for (n = 0; n<sz-1; n+=2) {
3224: i1 = idx[0];
3225: i2 = idx[1];
3226: idx += 2;
3227: tmp0 = x[i1];
3228: tmp1 = x[i2];
3229: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3230: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3231: }
3232: if (n == sz-1) {
3233: tmp0 = x[*idx];
3234: sum1 -= v1[0] * tmp0;
3235: sum2 -= v2[0] * tmp0;
3236: }
3237: x[row] = sum1*ibdiag[0] + sum2*ibdiag[2];
3238: x[row+1] = sum1*ibdiag[1] + sum2*ibdiag[3];
3239: break;
3240: case 3:
3241: v2 = a->a + ii[row+1];
3242: v3 = a->a + ii[row+2];
3243: sum1 = b[row];
3244: sum2 = b[row+1];
3245: sum3 = b[row+2];
3246: for (n = 0; n<sz-1; n+=2) {
3247: i1 = idx[0];
3248: i2 = idx[1];
3249: idx += 2;
3250: tmp0 = x[i1];
3251: tmp1 = x[i2];
3252: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3253: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3254: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3255: }
3256: if (n == sz-1) {
3257: tmp0 = x[*idx++];
3258: sum1 -= v1[0] * tmp0;
3259: sum2 -= v2[0] * tmp0;
3260: sum3 -= v3[0] * tmp0;
3261: v1++; v2++; v3++;
3262: }
3263: t[row] = sum1;
3264: t[row+1] = sum2;
3265: t[row+2] = sum3;
3266: sz = ii[row+1] - diag[row] - 3;
3267: idx = a->j + diag[row] + 3;
3268: v1 += 3;
3269: v2 += 3;
3270: v3 += 3;
3271: for (n = 0; n<sz-1; n+=2) {
3272: i1 = idx[0];
3273: i2 = idx[1];
3274: idx += 2;
3275: tmp0 = x[i1];
3276: tmp1 = x[i2];
3277: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3278: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3279: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3280: }
3281: if (n == sz-1) {
3282: tmp0 = x[*idx];
3283: sum1 -= v1[0] * tmp0;
3284: sum2 -= v2[0] * tmp0;
3285: sum3 -= v3[0] * tmp0;
3286: }
3287: x[row] = sum1*ibdiag[0] + sum2*ibdiag[3] + sum3*ibdiag[6];
3288: x[row+1] = sum1*ibdiag[1] + sum2*ibdiag[4] + sum3*ibdiag[7];
3289: x[row+2] = sum1*ibdiag[2] + sum2*ibdiag[5] + sum3*ibdiag[8];
3290: break;
3291: case 4:
3292: v2 = a->a + ii[row+1];
3293: v3 = a->a + ii[row+2];
3294: v4 = a->a + ii[row+3];
3295: sum1 = b[row];
3296: sum2 = b[row+1];
3297: sum3 = b[row+2];
3298: sum4 = b[row+3];
3299: for (n = 0; n<sz-1; n+=2) {
3300: i1 = idx[0];
3301: i2 = idx[1];
3302: idx += 2;
3303: tmp0 = x[i1];
3304: tmp1 = x[i2];
3305: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3306: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3307: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3308: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3309: }
3310: if (n == sz-1) {
3311: tmp0 = x[*idx++];
3312: sum1 -= v1[0] * tmp0;
3313: sum2 -= v2[0] * tmp0;
3314: sum3 -= v3[0] * tmp0;
3315: sum4 -= v4[0] * tmp0;
3316: v1++; v2++; v3++; v4++;
3317: }
3318: t[row] = sum1;
3319: t[row+1] = sum2;
3320: t[row+2] = sum3;
3321: t[row+3] = sum4;
3322: sz = ii[row+1] - diag[row] - 4;
3323: idx = a->j + diag[row] + 4;
3324: v1 += 4;
3325: v2 += 4;
3326: v3 += 4;
3327: v4 += 4;
3328: for (n = 0; n<sz-1; n+=2) {
3329: i1 = idx[0];
3330: i2 = idx[1];
3331: idx += 2;
3332: tmp0 = x[i1];
3333: tmp1 = x[i2];
3334: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3335: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3336: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3337: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3338: }
3339: if (n == sz-1) {
3340: tmp0 = x[*idx];
3341: sum1 -= v1[0] * tmp0;
3342: sum2 -= v2[0] * tmp0;
3343: sum3 -= v3[0] * tmp0;
3344: sum4 -= v4[0] * tmp0;
3345: }
3346: x[row] = sum1*ibdiag[0] + sum2*ibdiag[4] + sum3*ibdiag[8] + sum4*ibdiag[12];
3347: x[row+1] = sum1*ibdiag[1] + sum2*ibdiag[5] + sum3*ibdiag[9] + sum4*ibdiag[13];
3348: x[row+2] = sum1*ibdiag[2] + sum2*ibdiag[6] + sum3*ibdiag[10] + sum4*ibdiag[14];
3349: x[row+3] = sum1*ibdiag[3] + sum2*ibdiag[7] + sum3*ibdiag[11] + sum4*ibdiag[15];
3350: break;
3351: case 5:
3352: v2 = a->a + ii[row+1];
3353: v3 = a->a + ii[row+2];
3354: v4 = a->a + ii[row+3];
3355: v5 = a->a + ii[row+4];
3356: sum1 = b[row];
3357: sum2 = b[row+1];
3358: sum3 = b[row+2];
3359: sum4 = b[row+3];
3360: sum5 = b[row+4];
3361: for (n = 0; n<sz-1; n+=2) {
3362: i1 = idx[0];
3363: i2 = idx[1];
3364: idx += 2;
3365: tmp0 = x[i1];
3366: tmp1 = x[i2];
3367: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3368: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3369: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3370: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3371: sum5 -= v5[0] * tmp0 + v5[1] * tmp1; v5 += 2;
3372: }
3373: if (n == sz-1) {
3374: tmp0 = x[*idx++];
3375: sum1 -= v1[0] * tmp0;
3376: sum2 -= v2[0] * tmp0;
3377: sum3 -= v3[0] * tmp0;
3378: sum4 -= v4[0] * tmp0;
3379: sum5 -= v5[0] * tmp0;
3380: v1++; v2++; v3++; v4++; v5++;
3381: }
3382: t[row] = sum1;
3383: t[row+1] = sum2;
3384: t[row+2] = sum3;
3385: t[row+3] = sum4;
3386: t[row+4] = sum5;
3387: sz = ii[row+1] - diag[row] - 5;
3388: idx = a->j + diag[row] + 5;
3389: v1 += 5;
3390: v2 += 5;
3391: v3 += 5;
3392: v4 += 5;
3393: v5 += 5;
3394: for (n = 0; n<sz-1; n+=2) {
3395: i1 = idx[0];
3396: i2 = idx[1];
3397: idx += 2;
3398: tmp0 = x[i1];
3399: tmp1 = x[i2];
3400: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3401: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3402: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3403: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3404: sum5 -= v5[0] * tmp0 + v5[1] * tmp1; v5 += 2;
3405: }
3406: if (n == sz-1) {
3407: tmp0 = x[*idx];
3408: sum1 -= v1[0] * tmp0;
3409: sum2 -= v2[0] * tmp0;
3410: sum3 -= v3[0] * tmp0;
3411: sum4 -= v4[0] * tmp0;
3412: sum5 -= v5[0] * tmp0;
3413: }
3414: x[row] = sum1*ibdiag[0] + sum2*ibdiag[5] + sum3*ibdiag[10] + sum4*ibdiag[15] + sum5*ibdiag[20];
3415: x[row+1] = sum1*ibdiag[1] + sum2*ibdiag[6] + sum3*ibdiag[11] + sum4*ibdiag[16] + sum5*ibdiag[21];
3416: x[row+2] = sum1*ibdiag[2] + sum2*ibdiag[7] + sum3*ibdiag[12] + sum4*ibdiag[17] + sum5*ibdiag[22];
3417: x[row+3] = sum1*ibdiag[3] + sum2*ibdiag[8] + sum3*ibdiag[13] + sum4*ibdiag[18] + sum5*ibdiag[23];
3418: x[row+4] = sum1*ibdiag[4] + sum2*ibdiag[9] + sum3*ibdiag[14] + sum4*ibdiag[19] + sum5*ibdiag[24];
3419: break;
3420: default:
3421: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
3422: }
3423: }
3424: xb = t;
3425: PetscLogFlops(2.0*a->nz); /* undercounts diag inverse */
3426: } else xb = b;
3428: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
3430: ibdiag = a->inode.ibdiag+a->inode.bdiagsize;
3431: for (i=m-1, row=A->rmap->n-1; i>=0; i--) {
3432: ibdiag -= sizes[i]*sizes[i];
3434: /* set RHS */
3435: if (xb == b) {
3436: /* whole (old way) */
3437: sz = ii[row+1] - ii[row];
3438: idx = a->j + ii[row];
3439: switch (sizes[i]) {
3440: case 5:
3441: v5 = a->a + ii[row-4];
3442: case 4: /* fall through */
3443: v4 = a->a + ii[row-3];
3444: case 3:
3445: v3 = a->a + ii[row-2];
3446: case 2:
3447: v2 = a->a + ii[row-1];
3448: case 1:
3449: v1 = a->a + ii[row];
3450: break;
3451: default:
3452: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
3453: }
3454: } else {
3455: /* upper, no diag */
3456: sz = ii[row+1] - diag[row] - 1;
3457: idx = a->j + diag[row] + 1;
3458: switch (sizes[i]) {
3459: case 5:
3460: v5 = a->a + diag[row-4] + 5;
3461: case 4: /* fall through */
3462: v4 = a->a + diag[row-3] + 4;
3463: case 3:
3464: v3 = a->a + diag[row-2] + 3;
3465: case 2:
3466: v2 = a->a + diag[row-1] + 2;
3467: case 1:
3468: v1 = a->a + diag[row] + 1;
3469: }
3470: }
3471: /* set sum */
3472: switch (sizes[i]) {
3473: case 5:
3474: sum5 = xb[row-4];
3475: case 4: /* fall through */
3476: sum4 = xb[row-3];
3477: case 3:
3478: sum3 = xb[row-2];
3479: case 2:
3480: sum2 = xb[row-1];
3481: case 1:
3482: /* note that sum1 is associated with the last row */
3483: sum1 = xb[row];
3484: }
3485: /* do sums */
3486: for (n = 0; n<sz-1; n+=2) {
3487: i1 = idx[0];
3488: i2 = idx[1];
3489: idx += 2;
3490: tmp0 = x[i1];
3491: tmp1 = x[i2];
3492: switch (sizes[i]) {
3493: case 5:
3494: sum5 -= v5[0] * tmp0 + v5[1] * tmp1; v5 += 2;
3495: case 4: /* fall through */
3496: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3497: case 3:
3498: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3499: case 2:
3500: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3501: case 1:
3502: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3503: }
3504: }
3505: /* ragged edge */
3506: if (n == sz-1) {
3507: tmp0 = x[*idx];
3508: switch (sizes[i]) {
3509: case 5:
3510: sum5 -= *v5*tmp0;
3511: case 4: /* fall through */
3512: sum4 -= *v4*tmp0;
3513: case 3:
3514: sum3 -= *v3*tmp0;
3515: case 2:
3516: sum2 -= *v2*tmp0;
3517: case 1:
3518: sum1 -= *v1*tmp0;
3519: }
3520: }
3521: /* update */
3522: if (xb == b) {
3523: /* whole (old way) w/ diag */
3524: switch (sizes[i]) {
3525: case 5:
3526: x[row--] += sum5*ibdiag[4] + sum4*ibdiag[9] + sum3*ibdiag[14] + sum2*ibdiag[19] + sum1*ibdiag[24];
3527: x[row--] += sum5*ibdiag[3] + sum4*ibdiag[8] + sum3*ibdiag[13] + sum2*ibdiag[18] + sum1*ibdiag[23];
3528: x[row--] += sum5*ibdiag[2] + sum4*ibdiag[7] + sum3*ibdiag[12] + sum2*ibdiag[17] + sum1*ibdiag[22];
3529: x[row--] += sum5*ibdiag[1] + sum4*ibdiag[6] + sum3*ibdiag[11] + sum2*ibdiag[16] + sum1*ibdiag[21];
3530: x[row--] += sum5*ibdiag[0] + sum4*ibdiag[5] + sum3*ibdiag[10] + sum2*ibdiag[15] + sum1*ibdiag[20];
3531: break;
3532: case 4:
3533: x[row--] += sum4*ibdiag[3] + sum3*ibdiag[7] + sum2*ibdiag[11] + sum1*ibdiag[15];
3534: x[row--] += sum4*ibdiag[2] + sum3*ibdiag[6] + sum2*ibdiag[10] + sum1*ibdiag[14];
3535: x[row--] += sum4*ibdiag[1] + sum3*ibdiag[5] + sum2*ibdiag[9] + sum1*ibdiag[13];
3536: x[row--] += sum4*ibdiag[0] + sum3*ibdiag[4] + sum2*ibdiag[8] + sum1*ibdiag[12];
3537: break;
3538: case 3:
3539: x[row--] += sum3*ibdiag[2] + sum2*ibdiag[5] + sum1*ibdiag[8];
3540: x[row--] += sum3*ibdiag[1] + sum2*ibdiag[4] + sum1*ibdiag[7];
3541: x[row--] += sum3*ibdiag[0] + sum2*ibdiag[3] + sum1*ibdiag[6];
3542: break;
3543: case 2:
3544: x[row--] += sum2*ibdiag[1] + sum1*ibdiag[3];
3545: x[row--] += sum2*ibdiag[0] + sum1*ibdiag[2];
3546: break;
3547: case 1:
3548: x[row--] += sum1*(*ibdiag);
3549: break;
3550: }
3551: } else {
3552: /* no diag so set = */
3553: switch (sizes[i]) {
3554: case 5:
3555: x[row--] = sum5*ibdiag[4] + sum4*ibdiag[9] + sum3*ibdiag[14] + sum2*ibdiag[19] + sum1*ibdiag[24];
3556: x[row--] = sum5*ibdiag[3] + sum4*ibdiag[8] + sum3*ibdiag[13] + sum2*ibdiag[18] + sum1*ibdiag[23];
3557: x[row--] = sum5*ibdiag[2] + sum4*ibdiag[7] + sum3*ibdiag[12] + sum2*ibdiag[17] + sum1*ibdiag[22];
3558: x[row--] = sum5*ibdiag[1] + sum4*ibdiag[6] + sum3*ibdiag[11] + sum2*ibdiag[16] + sum1*ibdiag[21];
3559: x[row--] = sum5*ibdiag[0] + sum4*ibdiag[5] + sum3*ibdiag[10] + sum2*ibdiag[15] + sum1*ibdiag[20];
3560: break;
3561: case 4:
3562: x[row--] = sum4*ibdiag[3] + sum3*ibdiag[7] + sum2*ibdiag[11] + sum1*ibdiag[15];
3563: x[row--] = sum4*ibdiag[2] + sum3*ibdiag[6] + sum2*ibdiag[10] + sum1*ibdiag[14];
3564: x[row--] = sum4*ibdiag[1] + sum3*ibdiag[5] + sum2*ibdiag[9] + sum1*ibdiag[13];
3565: x[row--] = sum4*ibdiag[0] + sum3*ibdiag[4] + sum2*ibdiag[8] + sum1*ibdiag[12];
3566: break;
3567: case 3:
3568: x[row--] = sum3*ibdiag[2] + sum2*ibdiag[5] + sum1*ibdiag[8];
3569: x[row--] = sum3*ibdiag[1] + sum2*ibdiag[4] + sum1*ibdiag[7];
3570: x[row--] = sum3*ibdiag[0] + sum2*ibdiag[3] + sum1*ibdiag[6];
3571: break;
3572: case 2:
3573: x[row--] = sum2*ibdiag[1] + sum1*ibdiag[3];
3574: x[row--] = sum2*ibdiag[0] + sum1*ibdiag[2];
3575: break;
3576: case 1:
3577: x[row--] = sum1*(*ibdiag);
3578: break;
3579: }
3580: }
3581: }
3582: if (xb == b) {
3583: PetscLogFlops(2.0*a->nz);
3584: } else {
3585: PetscLogFlops(a->nz); /* assumes 1/2 in upper, undercounts diag inverse */
3586: }
3587: }
3588: }
3589: if (flag & SOR_EISENSTAT) {
3590: /*
3591: Apply (U + D)^-1 where D is now the block diagonal
3592: */
3593: ibdiag = a->inode.ibdiag+a->inode.bdiagsize;
3594: for (i=m-1, row=A->rmap->n-1; i>=0; i--) {
3595: ibdiag -= sizes[i]*sizes[i];
3596: sz = ii[row+1] - diag[row] - 1;
3597: v1 = a->a + diag[row] + 1;
3598: idx = a->j + diag[row] + 1;
3599: /* see comments for MatMult_SeqAIJ_Inode() for how this is coded */
3600: switch (sizes[i]) {
3601: case 1:
3603: sum1 = b[row];
3604: for (n = 0; n<sz-1; n+=2) {
3605: i1 = idx[0];
3606: i2 = idx[1];
3607: idx += 2;
3608: tmp0 = x[i1];
3609: tmp1 = x[i2];
3610: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3611: }
3613: if (n == sz-1) {
3614: tmp0 = x[*idx];
3615: sum1 -= *v1*tmp0;
3616: }
3617: x[row] = sum1*(*ibdiag);row--;
3618: break;
3620: case 2:
3622: sum1 = b[row];
3623: sum2 = b[row-1];
3624: /* note that sum1 is associated with the second of the two rows */
3625: v2 = a->a + diag[row-1] + 2;
3626: for (n = 0; n<sz-1; n+=2) {
3627: i1 = idx[0];
3628: i2 = idx[1];
3629: idx += 2;
3630: tmp0 = x[i1];
3631: tmp1 = x[i2];
3632: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3633: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3634: }
3636: if (n == sz-1) {
3637: tmp0 = x[*idx];
3638: sum1 -= *v1*tmp0;
3639: sum2 -= *v2*tmp0;
3640: }
3641: x[row] = sum2*ibdiag[1] + sum1*ibdiag[3];
3642: x[row-1] = sum2*ibdiag[0] + sum1*ibdiag[2];
3643: row -= 2;
3644: break;
3645: case 3:
3647: sum1 = b[row];
3648: sum2 = b[row-1];
3649: sum3 = b[row-2];
3650: v2 = a->a + diag[row-1] + 2;
3651: v3 = a->a + diag[row-2] + 3;
3652: for (n = 0; n<sz-1; n+=2) {
3653: i1 = idx[0];
3654: i2 = idx[1];
3655: idx += 2;
3656: tmp0 = x[i1];
3657: tmp1 = x[i2];
3658: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3659: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3660: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3661: }
3663: if (n == sz-1) {
3664: tmp0 = x[*idx];
3665: sum1 -= *v1*tmp0;
3666: sum2 -= *v2*tmp0;
3667: sum3 -= *v3*tmp0;
3668: }
3669: x[row] = sum3*ibdiag[2] + sum2*ibdiag[5] + sum1*ibdiag[8];
3670: x[row-1] = sum3*ibdiag[1] + sum2*ibdiag[4] + sum1*ibdiag[7];
3671: x[row-2] = sum3*ibdiag[0] + sum2*ibdiag[3] + sum1*ibdiag[6];
3672: row -= 3;
3673: break;
3674: case 4:
3676: sum1 = b[row];
3677: sum2 = b[row-1];
3678: sum3 = b[row-2];
3679: sum4 = b[row-3];
3680: v2 = a->a + diag[row-1] + 2;
3681: v3 = a->a + diag[row-2] + 3;
3682: v4 = a->a + diag[row-3] + 4;
3683: for (n = 0; n<sz-1; n+=2) {
3684: i1 = idx[0];
3685: i2 = idx[1];
3686: idx += 2;
3687: tmp0 = x[i1];
3688: tmp1 = x[i2];
3689: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3690: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3691: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3692: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3693: }
3695: if (n == sz-1) {
3696: tmp0 = x[*idx];
3697: sum1 -= *v1*tmp0;
3698: sum2 -= *v2*tmp0;
3699: sum3 -= *v3*tmp0;
3700: sum4 -= *v4*tmp0;
3701: }
3702: x[row] = sum4*ibdiag[3] + sum3*ibdiag[7] + sum2*ibdiag[11] + sum1*ibdiag[15];
3703: x[row-1] = sum4*ibdiag[2] + sum3*ibdiag[6] + sum2*ibdiag[10] + sum1*ibdiag[14];
3704: x[row-2] = sum4*ibdiag[1] + sum3*ibdiag[5] + sum2*ibdiag[9] + sum1*ibdiag[13];
3705: x[row-3] = sum4*ibdiag[0] + sum3*ibdiag[4] + sum2*ibdiag[8] + sum1*ibdiag[12];
3706: row -= 4;
3707: break;
3708: case 5:
3710: sum1 = b[row];
3711: sum2 = b[row-1];
3712: sum3 = b[row-2];
3713: sum4 = b[row-3];
3714: sum5 = b[row-4];
3715: v2 = a->a + diag[row-1] + 2;
3716: v3 = a->a + diag[row-2] + 3;
3717: v4 = a->a + diag[row-3] + 4;
3718: v5 = a->a + diag[row-4] + 5;
3719: for (n = 0; n<sz-1; n+=2) {
3720: i1 = idx[0];
3721: i2 = idx[1];
3722: idx += 2;
3723: tmp0 = x[i1];
3724: tmp1 = x[i2];
3725: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3726: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3727: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3728: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3729: sum5 -= v5[0] * tmp0 + v5[1] * tmp1; v5 += 2;
3730: }
3732: if (n == sz-1) {
3733: tmp0 = x[*idx];
3734: sum1 -= *v1*tmp0;
3735: sum2 -= *v2*tmp0;
3736: sum3 -= *v3*tmp0;
3737: sum4 -= *v4*tmp0;
3738: sum5 -= *v5*tmp0;
3739: }
3740: x[row] = sum5*ibdiag[4] + sum4*ibdiag[9] + sum3*ibdiag[14] + sum2*ibdiag[19] + sum1*ibdiag[24];
3741: x[row-1] = sum5*ibdiag[3] + sum4*ibdiag[8] + sum3*ibdiag[13] + sum2*ibdiag[18] + sum1*ibdiag[23];
3742: x[row-2] = sum5*ibdiag[2] + sum4*ibdiag[7] + sum3*ibdiag[12] + sum2*ibdiag[17] + sum1*ibdiag[22];
3743: x[row-3] = sum5*ibdiag[1] + sum4*ibdiag[6] + sum3*ibdiag[11] + sum2*ibdiag[16] + sum1*ibdiag[21];
3744: x[row-4] = sum5*ibdiag[0] + sum4*ibdiag[5] + sum3*ibdiag[10] + sum2*ibdiag[15] + sum1*ibdiag[20];
3745: row -= 5;
3746: break;
3747: default:
3748: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
3749: }
3750: }
3751: PetscLogFlops(a->nz);
3753: /*
3754: t = b - D x where D is the block diagonal
3755: */
3756: cnt = 0;
3757: for (i=0, row=0; i<m; i++) {
3758: switch (sizes[i]) {
3759: case 1:
3760: t[row] = b[row] - bdiag[cnt++]*x[row]; row++;
3761: break;
3762: case 2:
3763: x1 = x[row]; x2 = x[row+1];
3764: tmp1 = x1*bdiag[cnt] + x2*bdiag[cnt+2];
3765: tmp2 = x1*bdiag[cnt+1] + x2*bdiag[cnt+3];
3766: t[row] = b[row] - tmp1;
3767: t[row+1] = b[row+1] - tmp2; row += 2;
3768: cnt += 4;
3769: break;
3770: case 3:
3771: x1 = x[row]; x2 = x[row+1]; x3 = x[row+2];
3772: tmp1 = x1*bdiag[cnt] + x2*bdiag[cnt+3] + x3*bdiag[cnt+6];
3773: tmp2 = x1*bdiag[cnt+1] + x2*bdiag[cnt+4] + x3*bdiag[cnt+7];
3774: tmp3 = x1*bdiag[cnt+2] + x2*bdiag[cnt+5] + x3*bdiag[cnt+8];
3775: t[row] = b[row] - tmp1;
3776: t[row+1] = b[row+1] - tmp2;
3777: t[row+2] = b[row+2] - tmp3; row += 3;
3778: cnt += 9;
3779: break;
3780: case 4:
3781: x1 = x[row]; x2 = x[row+1]; x3 = x[row+2]; x4 = x[row+3];
3782: tmp1 = x1*bdiag[cnt] + x2*bdiag[cnt+4] + x3*bdiag[cnt+8] + x4*bdiag[cnt+12];
3783: tmp2 = x1*bdiag[cnt+1] + x2*bdiag[cnt+5] + x3*bdiag[cnt+9] + x4*bdiag[cnt+13];
3784: tmp3 = x1*bdiag[cnt+2] + x2*bdiag[cnt+6] + x3*bdiag[cnt+10] + x4*bdiag[cnt+14];
3785: tmp4 = x1*bdiag[cnt+3] + x2*bdiag[cnt+7] + x3*bdiag[cnt+11] + x4*bdiag[cnt+15];
3786: t[row] = b[row] - tmp1;
3787: t[row+1] = b[row+1] - tmp2;
3788: t[row+2] = b[row+2] - tmp3;
3789: t[row+3] = b[row+3] - tmp4; row += 4;
3790: cnt += 16;
3791: break;
3792: case 5:
3793: x1 = x[row]; x2 = x[row+1]; x3 = x[row+2]; x4 = x[row+3]; x5 = x[row+4];
3794: tmp1 = x1*bdiag[cnt] + x2*bdiag[cnt+5] + x3*bdiag[cnt+10] + x4*bdiag[cnt+15] + x5*bdiag[cnt+20];
3795: tmp2 = x1*bdiag[cnt+1] + x2*bdiag[cnt+6] + x3*bdiag[cnt+11] + x4*bdiag[cnt+16] + x5*bdiag[cnt+21];
3796: tmp3 = x1*bdiag[cnt+2] + x2*bdiag[cnt+7] + x3*bdiag[cnt+12] + x4*bdiag[cnt+17] + x5*bdiag[cnt+22];
3797: tmp4 = x1*bdiag[cnt+3] + x2*bdiag[cnt+8] + x3*bdiag[cnt+13] + x4*bdiag[cnt+18] + x5*bdiag[cnt+23];
3798: tmp5 = x1*bdiag[cnt+4] + x2*bdiag[cnt+9] + x3*bdiag[cnt+14] + x4*bdiag[cnt+19] + x5*bdiag[cnt+24];
3799: t[row] = b[row] - tmp1;
3800: t[row+1] = b[row+1] - tmp2;
3801: t[row+2] = b[row+2] - tmp3;
3802: t[row+3] = b[row+3] - tmp4;
3803: t[row+4] = b[row+4] - tmp5;row += 5;
3804: cnt += 25;
3805: break;
3806: default:
3807: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
3808: }
3809: }
3810: PetscLogFlops(m);
3814: /*
3815: Apply (L + D)^-1 where D is the block diagonal
3816: */
3817: for (i=0, row=0; i<m; i++) {
3818: sz = diag[row] - ii[row];
3819: v1 = a->a + ii[row];
3820: idx = a->j + ii[row];
3821: /* see comments for MatMult_SeqAIJ_Inode() for how this is coded */
3822: switch (sizes[i]) {
3823: case 1:
3825: sum1 = t[row];
3826: for (n = 0; n<sz-1; n+=2) {
3827: i1 = idx[0];
3828: i2 = idx[1];
3829: idx += 2;
3830: tmp0 = t[i1];
3831: tmp1 = t[i2];
3832: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3833: }
3835: if (n == sz-1) {
3836: tmp0 = t[*idx];
3837: sum1 -= *v1 * tmp0;
3838: }
3839: x[row] += t[row] = sum1*(*ibdiag++); row++;
3840: break;
3841: case 2:
3842: v2 = a->a + ii[row+1];
3843: sum1 = t[row];
3844: sum2 = t[row+1];
3845: for (n = 0; n<sz-1; n+=2) {
3846: i1 = idx[0];
3847: i2 = idx[1];
3848: idx += 2;
3849: tmp0 = t[i1];
3850: tmp1 = t[i2];
3851: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3852: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3853: }
3855: if (n == sz-1) {
3856: tmp0 = t[*idx];
3857: sum1 -= v1[0] * tmp0;
3858: sum2 -= v2[0] * tmp0;
3859: }
3860: x[row] += t[row] = sum1*ibdiag[0] + sum2*ibdiag[2];
3861: x[row+1] += t[row+1] = sum1*ibdiag[1] + sum2*ibdiag[3];
3862: ibdiag += 4; row += 2;
3863: break;
3864: case 3:
3865: v2 = a->a + ii[row+1];
3866: v3 = a->a + ii[row+2];
3867: sum1 = t[row];
3868: sum2 = t[row+1];
3869: sum3 = t[row+2];
3870: for (n = 0; n<sz-1; n+=2) {
3871: i1 = idx[0];
3872: i2 = idx[1];
3873: idx += 2;
3874: tmp0 = t[i1];
3875: tmp1 = t[i2];
3876: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3877: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3878: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3879: }
3881: if (n == sz-1) {
3882: tmp0 = t[*idx];
3883: sum1 -= v1[0] * tmp0;
3884: sum2 -= v2[0] * tmp0;
3885: sum3 -= v3[0] * tmp0;
3886: }
3887: x[row] += t[row] = sum1*ibdiag[0] + sum2*ibdiag[3] + sum3*ibdiag[6];
3888: x[row+1] += t[row+1] = sum1*ibdiag[1] + sum2*ibdiag[4] + sum3*ibdiag[7];
3889: x[row+2] += t[row+2] = sum1*ibdiag[2] + sum2*ibdiag[5] + sum3*ibdiag[8];
3890: ibdiag += 9; row += 3;
3891: break;
3892: case 4:
3893: v2 = a->a + ii[row+1];
3894: v3 = a->a + ii[row+2];
3895: v4 = a->a + ii[row+3];
3896: sum1 = t[row];
3897: sum2 = t[row+1];
3898: sum3 = t[row+2];
3899: sum4 = t[row+3];
3900: for (n = 0; n<sz-1; n+=2) {
3901: i1 = idx[0];
3902: i2 = idx[1];
3903: idx += 2;
3904: tmp0 = t[i1];
3905: tmp1 = t[i2];
3906: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3907: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3908: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3909: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3910: }
3912: if (n == sz-1) {
3913: tmp0 = t[*idx];
3914: sum1 -= v1[0] * tmp0;
3915: sum2 -= v2[0] * tmp0;
3916: sum3 -= v3[0] * tmp0;
3917: sum4 -= v4[0] * tmp0;
3918: }
3919: x[row] += t[row] = sum1*ibdiag[0] + sum2*ibdiag[4] + sum3*ibdiag[8] + sum4*ibdiag[12];
3920: x[row+1] += t[row+1] = sum1*ibdiag[1] + sum2*ibdiag[5] + sum3*ibdiag[9] + sum4*ibdiag[13];
3921: x[row+2] += t[row+2] = sum1*ibdiag[2] + sum2*ibdiag[6] + sum3*ibdiag[10] + sum4*ibdiag[14];
3922: x[row+3] += t[row+3] = sum1*ibdiag[3] + sum2*ibdiag[7] + sum3*ibdiag[11] + sum4*ibdiag[15];
3923: ibdiag += 16; row += 4;
3924: break;
3925: case 5:
3926: v2 = a->a + ii[row+1];
3927: v3 = a->a + ii[row+2];
3928: v4 = a->a + ii[row+3];
3929: v5 = a->a + ii[row+4];
3930: sum1 = t[row];
3931: sum2 = t[row+1];
3932: sum3 = t[row+2];
3933: sum4 = t[row+3];
3934: sum5 = t[row+4];
3935: for (n = 0; n<sz-1; n+=2) {
3936: i1 = idx[0];
3937: i2 = idx[1];
3938: idx += 2;
3939: tmp0 = t[i1];
3940: tmp1 = t[i2];
3941: sum1 -= v1[0] * tmp0 + v1[1] * tmp1; v1 += 2;
3942: sum2 -= v2[0] * tmp0 + v2[1] * tmp1; v2 += 2;
3943: sum3 -= v3[0] * tmp0 + v3[1] * tmp1; v3 += 2;
3944: sum4 -= v4[0] * tmp0 + v4[1] * tmp1; v4 += 2;
3945: sum5 -= v5[0] * tmp0 + v5[1] * tmp1; v5 += 2;
3946: }
3948: if (n == sz-1) {
3949: tmp0 = t[*idx];
3950: sum1 -= v1[0] * tmp0;
3951: sum2 -= v2[0] * tmp0;
3952: sum3 -= v3[0] * tmp0;
3953: sum4 -= v4[0] * tmp0;
3954: sum5 -= v5[0] * tmp0;
3955: }
3956: x[row] += t[row] = sum1*ibdiag[0] + sum2*ibdiag[5] + sum3*ibdiag[10] + sum4*ibdiag[15] + sum5*ibdiag[20];
3957: x[row+1] += t[row+1] = sum1*ibdiag[1] + sum2*ibdiag[6] + sum3*ibdiag[11] + sum4*ibdiag[16] + sum5*ibdiag[21];
3958: x[row+2] += t[row+2] = sum1*ibdiag[2] + sum2*ibdiag[7] + sum3*ibdiag[12] + sum4*ibdiag[17] + sum5*ibdiag[22];
3959: x[row+3] += t[row+3] = sum1*ibdiag[3] + sum2*ibdiag[8] + sum3*ibdiag[13] + sum4*ibdiag[18] + sum5*ibdiag[23];
3960: x[row+4] += t[row+4] = sum1*ibdiag[4] + sum2*ibdiag[9] + sum3*ibdiag[14] + sum4*ibdiag[19] + sum5*ibdiag[24];
3961: ibdiag += 25; row += 5;
3962: break;
3963: default:
3964: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
3965: }
3966: }
3967: PetscLogFlops(a->nz);
3968: }
3969: VecRestoreArray(xx,&x);
3970: VecRestoreArrayRead(bb,&b);
3971: return(0);
3972: }
3974: PetscErrorCode MatMultDiagonalBlock_SeqAIJ_Inode(Mat A,Vec bb,Vec xx)
3975: {
3976: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
3977: PetscScalar *x,tmp1,tmp2,tmp3,tmp4,tmp5,x1,x2,x3,x4,x5;
3978: const MatScalar *bdiag = a->inode.bdiag;
3979: const PetscScalar *b;
3980: PetscErrorCode ierr;
3981: PetscInt m = a->inode.node_count,cnt = 0,i,row;
3982: const PetscInt *sizes = a->inode.size;
3985: VecGetArray(xx,&x);
3986: VecGetArrayRead(bb,&b);
3987: cnt = 0;
3988: for (i=0, row=0; i<m; i++) {
3989: switch (sizes[i]) {
3990: case 1:
3991: x[row] = b[row]*bdiag[cnt++];row++;
3992: break;
3993: case 2:
3994: x1 = b[row]; x2 = b[row+1];
3995: tmp1 = x1*bdiag[cnt] + x2*bdiag[cnt+2];
3996: tmp2 = x1*bdiag[cnt+1] + x2*bdiag[cnt+3];
3997: x[row++] = tmp1;
3998: x[row++] = tmp2;
3999: cnt += 4;
4000: break;
4001: case 3:
4002: x1 = b[row]; x2 = b[row+1]; x3 = b[row+2];
4003: tmp1 = x1*bdiag[cnt] + x2*bdiag[cnt+3] + x3*bdiag[cnt+6];
4004: tmp2 = x1*bdiag[cnt+1] + x2*bdiag[cnt+4] + x3*bdiag[cnt+7];
4005: tmp3 = x1*bdiag[cnt+2] + x2*bdiag[cnt+5] + x3*bdiag[cnt+8];
4006: x[row++] = tmp1;
4007: x[row++] = tmp2;
4008: x[row++] = tmp3;
4009: cnt += 9;
4010: break;
4011: case 4:
4012: x1 = b[row]; x2 = b[row+1]; x3 = b[row+2]; x4 = b[row+3];
4013: tmp1 = x1*bdiag[cnt] + x2*bdiag[cnt+4] + x3*bdiag[cnt+8] + x4*bdiag[cnt+12];
4014: tmp2 = x1*bdiag[cnt+1] + x2*bdiag[cnt+5] + x3*bdiag[cnt+9] + x4*bdiag[cnt+13];
4015: tmp3 = x1*bdiag[cnt+2] + x2*bdiag[cnt+6] + x3*bdiag[cnt+10] + x4*bdiag[cnt+14];
4016: tmp4 = x1*bdiag[cnt+3] + x2*bdiag[cnt+7] + x3*bdiag[cnt+11] + x4*bdiag[cnt+15];
4017: x[row++] = tmp1;
4018: x[row++] = tmp2;
4019: x[row++] = tmp3;
4020: x[row++] = tmp4;
4021: cnt += 16;
4022: break;
4023: case 5:
4024: x1 = b[row]; x2 = b[row+1]; x3 = b[row+2]; x4 = b[row+3]; x5 = b[row+4];
4025: tmp1 = x1*bdiag[cnt] + x2*bdiag[cnt+5] + x3*bdiag[cnt+10] + x4*bdiag[cnt+15] + x5*bdiag[cnt+20];
4026: tmp2 = x1*bdiag[cnt+1] + x2*bdiag[cnt+6] + x3*bdiag[cnt+11] + x4*bdiag[cnt+16] + x5*bdiag[cnt+21];
4027: tmp3 = x1*bdiag[cnt+2] + x2*bdiag[cnt+7] + x3*bdiag[cnt+12] + x4*bdiag[cnt+17] + x5*bdiag[cnt+22];
4028: tmp4 = x1*bdiag[cnt+3] + x2*bdiag[cnt+8] + x3*bdiag[cnt+13] + x4*bdiag[cnt+18] + x5*bdiag[cnt+23];
4029: tmp5 = x1*bdiag[cnt+4] + x2*bdiag[cnt+9] + x3*bdiag[cnt+14] + x4*bdiag[cnt+19] + x5*bdiag[cnt+24];
4030: x[row++] = tmp1;
4031: x[row++] = tmp2;
4032: x[row++] = tmp3;
4033: x[row++] = tmp4;
4034: x[row++] = tmp5;
4035: cnt += 25;
4036: break;
4037: default:
4038: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Inode size %D not supported",sizes[i]);
4039: }
4040: }
4041: PetscLogFlops(2*cnt);
4042: VecRestoreArray(xx,&x);
4043: VecRestoreArrayRead(bb,&b);
4044: return(0);
4045: }
4047: /*
4048: samestructure indicates that the matrix has not changed its nonzero structure so we
4049: do not need to recompute the inodes
4050: */
4051: PetscErrorCode MatSeqAIJCheckInode(Mat A)
4052: {
4053: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4055: PetscInt i,j,m,nzx,nzy,*ns,node_count,blk_size;
4056: PetscBool flag;
4057: const PetscInt *idx,*idy,*ii;
4060: if (!a->inode.use) return(0);
4061: if (a->inode.checked && A->nonzerostate == a->inode.mat_nonzerostate) return(0);
4063: m = A->rmap->n;
4064: if (a->inode.size) ns = a->inode.size;
4065: else {
4066: PetscMalloc1(m+1,&ns);
4067: }
4069: i = 0;
4070: node_count = 0;
4071: idx = a->j;
4072: ii = a->i;
4073: while (i < m) { /* For each row */
4074: nzx = ii[i+1] - ii[i]; /* Number of nonzeros */
4075: /* Limits the number of elements in a node to 'a->inode.limit' */
4076: for (j=i+1,idy=idx,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) {
4077: nzy = ii[j+1] - ii[j]; /* Same number of nonzeros */
4078: if (nzy != nzx) break;
4079: idy += nzx; /* Same nonzero pattern */
4080: PetscMemcmp(idx,idy,nzx*sizeof(PetscInt),&flag);
4081: if (!flag) break;
4082: }
4083: ns[node_count++] = blk_size;
4084: idx += blk_size*nzx;
4085: i = j;
4086: }
4087: /* If not enough inodes found,, do not use inode version of the routines */
4088: if (!m || node_count > .8*m) {
4089: PetscFree(ns);
4091: a->inode.node_count = 0;
4092: a->inode.size = NULL;
4093: a->inode.use = PETSC_FALSE;
4094: A->ops->mult = MatMult_SeqAIJ;
4095: A->ops->sor = MatSOR_SeqAIJ;
4096: A->ops->multadd = MatMultAdd_SeqAIJ;
4097: A->ops->getrowij = MatGetRowIJ_SeqAIJ;
4098: A->ops->restorerowij = MatRestoreRowIJ_SeqAIJ;
4099: A->ops->getcolumnij = MatGetColumnIJ_SeqAIJ;
4100: A->ops->restorecolumnij = MatRestoreColumnIJ_SeqAIJ;
4101: A->ops->coloringpatch = 0;
4102: A->ops->multdiagonalblock = 0;
4104: PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);
4105: } else {
4106: if (!A->factortype) {
4107: A->ops->mult = MatMult_SeqAIJ_Inode;
4108: A->ops->sor = MatSOR_SeqAIJ_Inode;
4109: A->ops->multadd = MatMultAdd_SeqAIJ_Inode;
4110: A->ops->multdiagonalblock = MatMultDiagonalBlock_SeqAIJ_Inode;
4111: if (A->rmap->n == A->cmap->n) {
4112: A->ops->getrowij = MatGetRowIJ_SeqAIJ_Inode;
4113: A->ops->restorerowij = MatRestoreRowIJ_SeqAIJ_Inode;
4114: A->ops->getcolumnij = MatGetColumnIJ_SeqAIJ_Inode;
4115: A->ops->restorecolumnij = MatRestoreColumnIJ_SeqAIJ_Inode;
4116: A->ops->coloringpatch = MatColoringPatch_SeqAIJ_Inode;
4117: }
4118: } else {
4119: A->ops->solve = MatSolve_SeqAIJ_Inode_inplace;
4120: }
4121: a->inode.node_count = node_count;
4122: a->inode.size = ns;
4123: PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);
4124: }
4125: a->inode.checked = PETSC_TRUE;
4126: a->inode.mat_nonzerostate = A->nonzerostate;
4127: return(0);
4128: }
4130: PetscErrorCode MatDuplicate_SeqAIJ_Inode(Mat A,MatDuplicateOption cpvalues,Mat *C)
4131: {
4132: Mat B =*C;
4133: Mat_SeqAIJ *c=(Mat_SeqAIJ*)B->data,*a=(Mat_SeqAIJ*)A->data;
4135: PetscInt m=A->rmap->n;
4138: c->inode.use = a->inode.use;
4139: c->inode.limit = a->inode.limit;
4140: c->inode.max_limit = a->inode.max_limit;
4141: if (a->inode.size) {
4142: PetscMalloc1(m+1,&c->inode.size);
4143: c->inode.node_count = a->inode.node_count;
4144: PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(PetscInt));
4145: /* note the table of functions below should match that in MatSeqAIJCheckInode() */
4146: if (!B->factortype) {
4147: B->ops->mult = MatMult_SeqAIJ_Inode;
4148: B->ops->sor = MatSOR_SeqAIJ_Inode;
4149: B->ops->multadd = MatMultAdd_SeqAIJ_Inode;
4150: B->ops->getrowij = MatGetRowIJ_SeqAIJ_Inode;
4151: B->ops->restorerowij = MatRestoreRowIJ_SeqAIJ_Inode;
4152: B->ops->getcolumnij = MatGetColumnIJ_SeqAIJ_Inode;
4153: B->ops->restorecolumnij = MatRestoreColumnIJ_SeqAIJ_Inode;
4154: B->ops->coloringpatch = MatColoringPatch_SeqAIJ_Inode;
4155: B->ops->multdiagonalblock = MatMultDiagonalBlock_SeqAIJ_Inode;
4156: } else {
4157: B->ops->solve = MatSolve_SeqAIJ_Inode_inplace;
4158: }
4159: } else {
4160: c->inode.size = 0;
4161: c->inode.node_count = 0;
4162: }
4163: c->inode.ibdiagvalid = PETSC_FALSE;
4164: c->inode.ibdiag = 0;
4165: c->inode.bdiag = 0;
4166: return(0);
4167: }
4169: PETSC_STATIC_INLINE PetscErrorCode MatGetRow_FactoredLU(PetscInt *cols,PetscInt nzl,PetscInt nzu,PetscInt nz,const PetscInt *ai,const PetscInt *aj,const PetscInt *adiag,PetscInt row)
4170: {
4171: PetscInt k;
4172: const PetscInt *vi;
4175: vi = aj + ai[row];
4176: for (k=0; k<nzl; k++) cols[k] = vi[k];
4177: vi = aj + adiag[row];
4178: cols[nzl] = vi[0];
4179: vi = aj + adiag[row+1]+1;
4180: for (k=0; k<nzu; k++) cols[nzl+1+k] = vi[k];
4181: return(0);
4182: }
4183: /*
4184: MatSeqAIJCheckInode_FactorLU - Check Inode for factored seqaij matrix.
4185: Modified from MatSeqAIJCheckInode().
4187: Input Parameters:
4188: . Mat A - ILU or LU matrix factor
4190: */
4191: PetscErrorCode MatSeqAIJCheckInode_FactorLU(Mat A)
4192: {
4193: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4195: PetscInt i,j,m,nzl1,nzu1,nzl2,nzu2,nzx,nzy,node_count,blk_size;
4196: PetscInt *cols1,*cols2,*ns;
4197: const PetscInt *ai = a->i,*aj = a->j, *adiag = a->diag;
4198: PetscBool flag;
4201: if (!a->inode.use) return(0);
4202: if (a->inode.checked) return(0);
4204: m = A->rmap->n;
4205: if (a->inode.size) ns = a->inode.size;
4206: else {
4207: PetscMalloc1(m+1,&ns);
4208: }
4210: i = 0;
4211: node_count = 0;
4212: PetscMalloc2(m,&cols1,m,&cols2);
4213: while (i < m) { /* For each row */
4214: nzl1 = ai[i+1] - ai[i]; /* Number of nonzeros in L */
4215: nzu1 = adiag[i] - adiag[i+1] - 1; /* Number of nonzeros in U excluding diagonal*/
4216: nzx = nzl1 + nzu1 + 1;
4217: MatGetRow_FactoredLU(cols1,nzl1,nzu1,nzx,ai,aj,adiag,i);
4219: /* Limits the number of elements in a node to 'a->inode.limit' */
4220: for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) {
4221: nzl2 = ai[j+1] - ai[j];
4222: nzu2 = adiag[j] - adiag[j+1] - 1;
4223: nzy = nzl2 + nzu2 + 1;
4224: if (nzy != nzx) break;
4225: MatGetRow_FactoredLU(cols2,nzl2,nzu2,nzy,ai,aj,adiag,j);
4226: PetscMemcmp(cols1,cols2,nzx*sizeof(PetscInt),&flag);
4227: if (!flag) break;
4228: }
4229: ns[node_count++] = blk_size;
4230: i = j;
4231: }
4232: PetscFree2(cols1,cols2);
4233: /* If not enough inodes found,, do not use inode version of the routines */
4234: if (!m || node_count > .8*m) {
4235: PetscFree(ns);
4237: a->inode.node_count = 0;
4238: a->inode.size = NULL;
4239: a->inode.use = PETSC_FALSE;
4241: PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);
4242: } else {
4243: A->ops->mult = 0;
4244: A->ops->sor = 0;
4245: A->ops->multadd = 0;
4246: A->ops->getrowij = 0;
4247: A->ops->restorerowij = 0;
4248: A->ops->getcolumnij = 0;
4249: A->ops->restorecolumnij = 0;
4250: A->ops->coloringpatch = 0;
4251: A->ops->multdiagonalblock = 0;
4252: a->inode.node_count = node_count;
4253: a->inode.size = ns;
4255: PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);
4256: }
4257: a->inode.checked = PETSC_TRUE;
4258: return(0);
4259: }
4261: PetscErrorCode MatSeqAIJInvalidateDiagonal_Inode(Mat A)
4262: {
4263: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
4266: a->inode.ibdiagvalid = PETSC_FALSE;
4267: return(0);
4268: }
4270: /*
4271: This is really ugly. if inodes are used this replaces the
4272: permutations with ones that correspond to rows/cols of the matrix
4273: rather then inode blocks
4274: */
4275: PetscErrorCode MatInodeAdjustForInodes(Mat A,IS *rperm,IS *cperm)
4276: {
4280: PetscTryMethod(A,"MatInodeAdjustForInodes_C",(Mat,IS*,IS*),(A,rperm,cperm));
4281: return(0);
4282: }
4284: PetscErrorCode MatInodeAdjustForInodes_SeqAIJ_Inode(Mat A,IS *rperm,IS *cperm)
4285: {
4286: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
4288: PetscInt m = A->rmap->n,n = A->cmap->n,i,j,nslim_row = a->inode.node_count;
4289: const PetscInt *ridx,*cidx;
4290: PetscInt row,col,*permr,*permc,*ns_row = a->inode.size,*tns,start_val,end_val,indx;
4291: PetscInt nslim_col,*ns_col;
4292: IS ris = *rperm,cis = *cperm;
4295: if (!a->inode.size) return(0); /* no inodes so return */
4296: if (a->inode.node_count == m) return(0); /* all inodes are of size 1 */
4298: MatCreateColInode_Private(A,&nslim_col,&ns_col);
4299: PetscMalloc1(((nslim_row>nslim_col) ? nslim_row : nslim_col)+1,&tns);
4300: PetscMalloc2(m,&permr,n,&permc);
4302: ISGetIndices(ris,&ridx);
4303: ISGetIndices(cis,&cidx);
4305: /* Form the inode structure for the rows of permuted matric using inv perm*/
4306: for (i=0,tns[0]=0; i<nslim_row; ++i) tns[i+1] = tns[i] + ns_row[i];
4308: /* Construct the permutations for rows*/
4309: for (i=0,row = 0; i<nslim_row; ++i) {
4310: indx = ridx[i];
4311: start_val = tns[indx];
4312: end_val = tns[indx + 1];
4313: for (j=start_val; j<end_val; ++j,++row) permr[row]= j;
4314: }
4316: /* Form the inode structure for the columns of permuted matrix using inv perm*/
4317: for (i=0,tns[0]=0; i<nslim_col; ++i) tns[i+1] = tns[i] + ns_col[i];
4319: /* Construct permutations for columns */
4320: for (i=0,col=0; i<nslim_col; ++i) {
4321: indx = cidx[i];
4322: start_val = tns[indx];
4323: end_val = tns[indx + 1];
4324: for (j = start_val; j<end_val; ++j,++col) permc[col]= j;
4325: }
4327: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,rperm);
4328: ISSetPermutation(*rperm);
4329: ISCreateGeneral(PETSC_COMM_SELF,n,permc,PETSC_COPY_VALUES,cperm);
4330: ISSetPermutation(*cperm);
4332: ISRestoreIndices(ris,&ridx);
4333: ISRestoreIndices(cis,&cidx);
4335: PetscFree(ns_col);
4336: PetscFree2(permr,permc);
4337: ISDestroy(&cis);
4338: ISDestroy(&ris);
4339: PetscFree(tns);
4340: return(0);
4341: }
4343: /*@C
4344: MatInodeGetInodeSizes - Returns the inode information of the Inode matrix.
4346: Not Collective
4348: Input Parameter:
4349: . A - the Inode matrix or matrix derived from the Inode class -- e.g., SeqAIJ
4351: Output Parameter:
4352: + node_count - no of inodes present in the matrix.
4353: . sizes - an array of size node_count,with sizes of each inode.
4354: - limit - the max size used to generate the inodes.
4356: Level: advanced
4358: Notes:
4359: This routine returns some internal storage information
4360: of the matrix, it is intended to be used by advanced users.
4361: It should be called after the matrix is assembled.
4362: The contents of the sizes[] array should not be changed.
4363: NULL may be passed for information not requested.
4365: .keywords: matrix, seqaij, get, inode
4367: .seealso: MatGetInfo()
4368: @*/
4369: PetscErrorCode MatInodeGetInodeSizes(Mat A,PetscInt *node_count,PetscInt *sizes[],PetscInt *limit)
4370: {
4371: PetscErrorCode ierr,(*f)(Mat,PetscInt*,PetscInt*[],PetscInt*);
4374: if (!A->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4375: PetscObjectQueryFunction((PetscObject)A,"MatInodeGetInodeSizes_C",&f);
4376: if (f) {
4377: (*f)(A,node_count,sizes,limit);
4378: }
4379: return(0);
4380: }
4382: PetscErrorCode MatInodeGetInodeSizes_SeqAIJ_Inode(Mat A,PetscInt *node_count,PetscInt *sizes[],PetscInt *limit)
4383: {
4384: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
4387: if (node_count) *node_count = a->inode.node_count;
4388: if (sizes) *sizes = a->inode.size;
4389: if (limit) *limit = a->inode.limit;
4390: return(0);
4391: }