Actual source code: baijsolvnat7.c
petsc-3.12.5 2020-03-29
1: #include <../src/mat/impls/baij/seq/baij.h>
2: #include <petsc/private/kernels/blockinvert.h>
4: PetscErrorCode MatSolve_SeqBAIJ_7_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx)
5: {
6: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
7: const PetscInt *diag=a->diag,n=a->mbs,*vi,*ai=a->i,*aj=a->j;
8: PetscErrorCode ierr;
9: PetscInt i,nz,idx,idt,jdx;
10: const MatScalar *aa=a->a,*v;
11: PetscScalar *x,s1,s2,s3,s4,s5,s6,s7,x1,x2,x3,x4,x5,x6,x7;
12: const PetscScalar *b;
15: VecGetArrayRead(bb,&b);
16: VecGetArray(xx,&x);
17: /* forward solve the lower triangular */
18: idx = 0;
19: x[0] = b[idx]; x[1] = b[1+idx]; x[2] = b[2+idx];
20: x[3] = b[3+idx]; x[4] = b[4+idx]; x[5] = b[5+idx];
21: x[6] = b[6+idx];
22: for (i=1; i<n; i++) {
23: v = aa + 49*ai[i];
24: vi = aj + ai[i];
25: nz = diag[i] - ai[i];
26: idx = 7*i;
27: s1 = b[idx]; s2 = b[1+idx]; s3 = b[2+idx];
28: s4 = b[3+idx]; s5 = b[4+idx]; s6 = b[5+idx];
29: s7 = b[6+idx];
30: while (nz--) {
31: jdx = 7*(*vi++);
32: x1 = x[jdx]; x2 = x[1+jdx]; x3 = x[2+jdx];
33: x4 = x[3+jdx]; x5 = x[4+jdx]; x6 = x[5+jdx];
34: x7 = x[6+jdx];
35: s1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7;
36: s2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7;
37: s3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7;
38: s4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7;
39: s5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7;
40: s6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7;
41: s7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7;
42: v += 49;
43: }
44: x[idx] = s1;
45: x[1+idx] = s2;
46: x[2+idx] = s3;
47: x[3+idx] = s4;
48: x[4+idx] = s5;
49: x[5+idx] = s6;
50: x[6+idx] = s7;
51: }
52: /* backward solve the upper triangular */
53: for (i=n-1; i>=0; i--) {
54: v = aa + 49*diag[i] + 49;
55: vi = aj + diag[i] + 1;
56: nz = ai[i+1] - diag[i] - 1;
57: idt = 7*i;
58: s1 = x[idt]; s2 = x[1+idt];
59: s3 = x[2+idt]; s4 = x[3+idt];
60: s5 = x[4+idt]; s6 = x[5+idt];
61: s7 = x[6+idt];
62: while (nz--) {
63: idx = 7*(*vi++);
64: x1 = x[idx]; x2 = x[1+idx]; x3 = x[2+idx];
65: x4 = x[3+idx]; x5 = x[4+idx]; x6 = x[5+idx];
66: x7 = x[6+idx];
67: s1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7;
68: s2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7;
69: s3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7;
70: s4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7;
71: s5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7;
72: s6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7;
73: s7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7;
74: v += 49;
75: }
76: v = aa + 49*diag[i];
77: x[idt] = v[0]*s1 + v[7]*s2 + v[14]*s3 + v[21]*s4
78: + v[28]*s5 + v[35]*s6 + v[42]*s7;
79: x[1+idt] = v[1]*s1 + v[8]*s2 + v[15]*s3 + v[22]*s4
80: + v[29]*s5 + v[36]*s6 + v[43]*s7;
81: x[2+idt] = v[2]*s1 + v[9]*s2 + v[16]*s3 + v[23]*s4
82: + v[30]*s5 + v[37]*s6 + v[44]*s7;
83: x[3+idt] = v[3]*s1 + v[10]*s2 + v[17]*s3 + v[24]*s4
84: + v[31]*s5 + v[38]*s6 + v[45]*s7;
85: x[4+idt] = v[4]*s1 + v[11]*s2 + v[18]*s3 + v[25]*s4
86: + v[32]*s5 + v[39]*s6 + v[46]*s7;
87: x[5+idt] = v[5]*s1 + v[12]*s2 + v[19]*s3 + v[26]*s4
88: + v[33]*s5 + v[40]*s6 + v[47]*s7;
89: x[6+idt] = v[6]*s1 + v[13]*s2 + v[20]*s3 + v[27]*s4
90: + v[34]*s5 + v[41]*s6 + v[48]*s7;
91: }
93: VecRestoreArrayRead(bb,&b);
94: VecRestoreArray(xx,&x);
95: PetscLogFlops(2.0*36*(a->nz) - 6.0*A->cmap->n);
96: return(0);
97: }
99: PetscErrorCode MatSolve_SeqBAIJ_7_NaturalOrdering(Mat A,Vec bb,Vec xx)
100: {
101: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
102: const PetscInt n =a->mbs,*vi,*ai=a->i,*aj=a->j,*adiag=a->diag;
103: PetscErrorCode ierr;
104: PetscInt i,k,nz,idx,jdx,idt;
105: const PetscInt bs = A->rmap->bs,bs2 = a->bs2;
106: const MatScalar *aa=a->a,*v;
107: PetscScalar *x;
108: const PetscScalar *b;
109: PetscScalar s1,s2,s3,s4,s5,s6,s7,x1,x2,x3,x4,x5,x6,x7;
112: VecGetArrayRead(bb,&b);
113: VecGetArray(xx,&x);
114: /* forward solve the lower triangular */
115: idx = 0;
116: x[0] = b[idx]; x[1] = b[1+idx];x[2] = b[2+idx];x[3] = b[3+idx];
117: x[4] = b[4+idx];x[5] = b[5+idx];x[6] = b[6+idx];
118: for (i=1; i<n; i++) {
119: v = aa + bs2*ai[i];
120: vi = aj + ai[i];
121: nz = ai[i+1] - ai[i];
122: idx = bs*i;
123: s1 = b[idx];s2 = b[1+idx];s3 = b[2+idx];s4 = b[3+idx];
124: s5 = b[4+idx];s6 = b[5+idx];s7 = b[6+idx];
125: for (k=0; k<nz; k++) {
126: jdx = bs*vi[k];
127: x1 = x[jdx];x2 = x[1+jdx]; x3 =x[2+jdx];x4 =x[3+jdx];
128: x5 = x[4+jdx]; x6 = x[5+jdx];x7 = x[6+jdx];
129: s1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7;
130: s2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7;
131: s3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7;
132: s4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7;
133: s5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7;
134: s6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7;
135: s7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7;
136: v += bs2;
137: }
139: x[idx] = s1;
140: x[1+idx] = s2;
141: x[2+idx] = s3;
142: x[3+idx] = s4;
143: x[4+idx] = s5;
144: x[5+idx] = s6;
145: x[6+idx] = s7;
146: }
148: /* backward solve the upper triangular */
149: for (i=n-1; i>=0; i--) {
150: v = aa + bs2*(adiag[i+1]+1);
151: vi = aj + adiag[i+1]+1;
152: nz = adiag[i] - adiag[i+1]-1;
153: idt = bs*i;
154: s1 = x[idt]; s2 = x[1+idt];s3 = x[2+idt];s4 = x[3+idt];
155: s5 = x[4+idt];s6 = x[5+idt];s7 = x[6+idt];
156: for (k=0; k<nz; k++) {
157: idx = bs*vi[k];
158: x1 = x[idx]; x2 = x[1+idx]; x3 = x[2+idx];x4 = x[3+idx];
159: x5 = x[4+idx];x6 = x[5+idx];x7 = x[6+idx];
160: s1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7;
161: s2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7;
162: s3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7;
163: s4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7;
164: s5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7;
165: s6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7;
166: s7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7;
167: v += bs2;
168: }
169: /* x = inv_diagonal*x */
170: x[idt] = v[0]*s1 + v[7]*s2 + v[14]*s3 + v[21]*s4 + v[28]*s5 + v[35]*s6 + v[42]*s7;
171: x[1+idt] = v[1]*s1 + v[8]*s2 + v[15]*s3 + v[22]*s4 + v[29]*s5 + v[36]*s6 + v[43]*s7;
172: x[2+idt] = v[2]*s1 + v[9]*s2 + v[16]*s3 + v[23]*s4 + v[30]*s5 + v[37]*s6 + v[44]*s7;
173: x[3+idt] = v[3]*s1 + v[10]*s2 + v[17]*s3 + v[24]*s4 + v[31]*s5 + v[38]*s6 + v[45]*s7;
174: x[4+idt] = v[4]*s1 + v[11]*s2 + v[18]*s3 + v[25]*s4 + v[32]*s5 + v[39]*s6 + v[46]*s7;
175: x[5+idt] = v[5]*s1 + v[12]*s2 + v[19]*s3 + v[26]*s4 + v[33]*s5 + v[40]*s6 + v[47]*s7;
176: x[6+idt] = v[6]*s1 + v[13]*s2 + v[20]*s3 + v[27]*s4 + v[34]*s5 + v[41]*s6 + v[48]*s7;
177: }
179: VecRestoreArrayRead(bb,&b);
180: VecRestoreArray(xx,&x);
181: PetscLogFlops(2.0*bs2*(a->nz) - bs*A->cmap->n);
182: return(0);
183: }