Actual source code: baijsolvnat7.c

petsc-master 2019-06-22
Report Typos and Errors
  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: }