Actual source code: baijsolvnat6.c

  1: #include <../src/mat/impls/baij/seq/baij.h>
  2: #include <petsc/private/kernels/blockinvert.h>

  4: PetscErrorCode MatSolve_SeqBAIJ_6_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx)
  5: {
  6:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
  7:   PetscInt          i,nz,idx,idt,jdx;
  8:   const PetscInt    *diag = a->diag,*vi,n=a->mbs,*ai=a->i,*aj=a->j;
  9:   const MatScalar   *aa   =a->a,*v;
 10:   PetscScalar       *x,s1,s2,s3,s4,s5,s6,x1,x2,x3,x4,x5,x6;
 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:   for (i=1; i<n; i++) {
 20:     v   =  aa + 36*ai[i];
 21:     vi  =  aj + ai[i];
 22:     nz  =  diag[i] - ai[i];
 23:     idx =  6*i;
 24:     s1  =  b[idx];   s2 = b[1+idx]; s3 = b[2+idx];
 25:     s4  =  b[3+idx]; s5 = b[4+idx]; s6 = b[5+idx];
 26:     while (nz--) {
 27:       jdx = 6*(*vi++);
 28:       x1  = x[jdx];   x2 = x[1+jdx]; x3 = x[2+jdx];
 29:       x4  = x[3+jdx]; x5 = x[4+jdx]; x6 = x[5+jdx];
 30:       s1 -= v[0]*x1 + v[6]*x2  + v[12]*x3 + v[18]*x4 + v[24]*x5 + v[30]*x6;
 31:       s2 -= v[1]*x1 + v[7]*x2  + v[13]*x3 + v[19]*x4 + v[25]*x5 + v[31]*x6;
 32:       s3 -= v[2]*x1 + v[8]*x2  + v[14]*x3 + v[20]*x4 + v[26]*x5 + v[32]*x6;
 33:       s4 -= v[3]*x1 + v[9]*x2  + v[15]*x3 + v[21]*x4 + v[27]*x5 + v[33]*x6;
 34:       s5 -= v[4]*x1 + v[10]*x2 + v[16]*x3 + v[22]*x4 + v[28]*x5 + v[34]*x6;
 35:       s6 -= v[5]*x1 + v[11]*x2 + v[17]*x3 + v[23]*x4 + v[29]*x5 + v[35]*x6;
 36:       v  += 36;
 37:     }
 38:     x[idx]   = s1;
 39:     x[1+idx] = s2;
 40:     x[2+idx] = s3;
 41:     x[3+idx] = s4;
 42:     x[4+idx] = s5;
 43:     x[5+idx] = s6;
 44:   }
 45:   /* backward solve the upper triangular */
 46:   for (i=n-1; i>=0; i--) {
 47:     v   = aa + 36*diag[i] + 36;
 48:     vi  = aj + diag[i] + 1;
 49:     nz  = ai[i+1] - diag[i] - 1;
 50:     idt = 6*i;
 51:     s1  = x[idt];   s2 = x[1+idt];
 52:     s3  = x[2+idt]; s4 = x[3+idt];
 53:     s5  = x[4+idt]; s6 = x[5+idt];
 54:     while (nz--) {
 55:       idx = 6*(*vi++);
 56:       x1  = x[idx];   x2 = x[1+idx]; x3 = x[2+idx];
 57:       x4  = x[3+idx]; x5 = x[4+idx]; x6 = x[5+idx];
 58:       s1 -= v[0]*x1 + v[6]*x2  + v[12]*x3 + v[18]*x4 + v[24]*x5 + v[30]*x6;
 59:       s2 -= v[1]*x1 + v[7]*x2  + v[13]*x3 + v[19]*x4 + v[25]*x5 + v[31]*x6;
 60:       s3 -= v[2]*x1 + v[8]*x2  + v[14]*x3 + v[20]*x4 + v[26]*x5 + v[32]*x6;
 61:       s4 -= v[3]*x1 + v[9]*x2  + v[15]*x3 + v[21]*x4 + v[27]*x5 + v[33]*x6;
 62:       s5 -= v[4]*x1 + v[10]*x2 + v[16]*x3 + v[22]*x4 + v[28]*x5 + v[34]*x6;
 63:       s6 -= v[5]*x1 + v[11]*x2 + v[17]*x3 + v[23]*x4 + v[29]*x5 + v[35]*x6;
 64:       v  += 36;
 65:     }
 66:     v        = aa + 36*diag[i];
 67:     x[idt]   = v[0]*s1 + v[6]*s2  + v[12]*s3 + v[18]*s4 + v[24]*s5 + v[30]*s6;
 68:     x[1+idt] = v[1]*s1 + v[7]*s2  + v[13]*s3 + v[19]*s4 + v[25]*s5 + v[31]*s6;
 69:     x[2+idt] = v[2]*s1 + v[8]*s2  + v[14]*s3 + v[20]*s4 + v[26]*s5 + v[32]*s6;
 70:     x[3+idt] = v[3]*s1 + v[9]*s2  + v[15]*s3 + v[21]*s4 + v[27]*s5 + v[33]*s6;
 71:     x[4+idt] = v[4]*s1 + v[10]*s2 + v[16]*s3 + v[22]*s4 + v[28]*s5 + v[34]*s6;
 72:     x[5+idt] = v[5]*s1 + v[11]*s2 + v[17]*s3 + v[23]*s4 + v[29]*s5 + v[35]*s6;
 73:   }

 75:   VecRestoreArrayRead(bb,&b);
 76:   VecRestoreArray(xx,&x);
 77:   PetscLogFlops(2.0*36*(a->nz) - 6.0*A->cmap->n);
 78:   return 0;
 79: }

 81: PetscErrorCode MatSolve_SeqBAIJ_6_NaturalOrdering(Mat A,Vec bb,Vec xx)
 82: {
 83:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
 84:   const PetscInt    n  =a->mbs,*vi,*ai=a->i,*aj=a->j,*adiag=a->diag;
 85:   PetscInt          i,k,nz,idx,jdx,idt;
 86:   const PetscInt    bs = A->rmap->bs,bs2 = a->bs2;
 87:   const MatScalar   *aa=a->a,*v;
 88:   PetscScalar       *x;
 89:   const PetscScalar *b;
 90:   PetscScalar       s1,s2,s3,s4,s5,s6,x1,x2,x3,x4,x5,x6;

 92:   VecGetArrayRead(bb,&b);
 93:   VecGetArray(xx,&x);
 94:   /* forward solve the lower triangular */
 95:   idx  = 0;
 96:   x[0] = b[idx]; x[1] = b[1+idx];x[2] = b[2+idx];x[3] = b[3+idx];
 97:   x[4] = b[4+idx];x[5] = b[5+idx];
 98:   for (i=1; i<n; i++) {
 99:     v   = aa + bs2*ai[i];
100:     vi  = aj + ai[i];
101:     nz  = ai[i+1] - ai[i];
102:     idx = bs*i;
103:     s1  = b[idx];s2 = b[1+idx];s3 = b[2+idx];s4 = b[3+idx];
104:     s5  = b[4+idx];s6 = b[5+idx];
105:     for (k=0; k<nz; k++) {
106:       jdx = bs*vi[k];
107:       x1  = x[jdx];x2 = x[1+jdx]; x3 =x[2+jdx];x4 =x[3+jdx];
108:       x5  = x[4+jdx]; x6 = x[5+jdx];
109:       s1 -= v[0]*x1 + v[6]*x2 + v[12]*x3 + v[18]*x4  + v[24]*x5 + v[30]*x6;
110:       s2 -= v[1]*x1 + v[7]*x2 + v[13]*x3 + v[19]*x4  + v[25]*x5 + v[31]*x6;
111:       s3 -= v[2]*x1 + v[8]*x2 + v[14]*x3 + v[20]*x4  + v[26]*x5 + v[32]*x6;
112:       s4 -= v[3]*x1 + v[9]*x2 + v[15]*x3 + v[21]*x4  + v[27]*x5 + v[33]*x6;
113:       s5 -= v[4]*x1 + v[10]*x2 + v[16]*x3 + v[22]*x4  + v[28]*x5 + v[34]*x6;
114:       s6 -= v[5]*x1 + v[11]*x2 + v[17]*x3 + v[23]*x4  + v[29]*x5 + v[35]*x6;
115:       v  +=  bs2;
116:     }

118:     x[idx]   = s1;
119:     x[1+idx] = s2;
120:     x[2+idx] = s3;
121:     x[3+idx] = s4;
122:     x[4+idx] = s5;
123:     x[5+idx] = s6;
124:   }

126:   /* backward solve the upper triangular */
127:   for (i=n-1; i>=0; i--) {
128:     v   = aa + bs2*(adiag[i+1]+1);
129:     vi  = aj + adiag[i+1]+1;
130:     nz  = adiag[i] - adiag[i+1]-1;
131:     idt = bs*i;
132:     s1  = x[idt];  s2 = x[1+idt];s3 = x[2+idt];s4 = x[3+idt];
133:     s5  = x[4+idt];s6 = x[5+idt];
134:     for (k=0; k<nz; k++) {
135:       idx = bs*vi[k];
136:       x1  = x[idx];   x2 = x[1+idx]; x3 = x[2+idx];x4 = x[3+idx];
137:       x5  = x[4+idx];x6 = x[5+idx];
138:       s1 -= v[0]*x1 + v[6]*x2 + v[12]*x3 + v[18]*x4  + v[24]*x5 + v[30]*x6;
139:       s2 -= v[1]*x1 + v[7]*x2 + v[13]*x3 + v[19]*x4  + v[25]*x5 + v[31]*x6;
140:       s3 -= v[2]*x1 + v[8]*x2 + v[14]*x3 + v[20]*x4  + v[26]*x5 + v[32]*x6;
141:       s4 -= v[3]*x1 + v[9]*x2 + v[15]*x3 + v[21]*x4  + v[27]*x5 + v[33]*x6;
142:       s5 -= v[4]*x1 + v[10]*x2 + v[16]*x3 + v[22]*x4  + v[28]*x5 + v[34]*x6;
143:       s6 -= v[5]*x1 + v[11]*x2 + v[17]*x3 + v[23]*x4  + v[29]*x5 + v[35]*x6;
144:       v  +=  bs2;
145:     }
146:     /* x = inv_diagonal*x */
147:     x[idt]   = v[0]*s1 + v[6]*s2 + v[12]*s3 + v[18]*s4 + v[24]*s5 + v[30]*s6;
148:     x[1+idt] = v[1]*s1 + v[7]*s2 + v[13]*s3 + v[19]*s4 + v[25]*s5 + v[31]*s6;
149:     x[2+idt] = v[2]*s1 + v[8]*s2 + v[14]*s3 + v[20]*s4 + v[26]*s5 + v[32]*s6;
150:     x[3+idt] = v[3]*s1 + v[9]*s2 + v[15]*s3 + v[21]*s4 + v[27]*s5 + v[33]*s6;
151:     x[4+idt] = v[4]*s1 + v[10]*s2 + v[16]*s3 + v[22]*s4 + v[28]*s5 + v[34]*s6;
152:     x[5+idt] = v[5]*s1 + v[11]*s2 + v[17]*s3 + v[23]*s4 + v[29]*s5 + v[35]*s6;
153:   }

155:   VecRestoreArrayRead(bb,&b);
156:   VecRestoreArray(xx,&x);
157:   PetscLogFlops(2.0*bs2*(a->nz) - bs*A->cmap->n);
158:   return 0;
159: }