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