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