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