Actual source code: baijsolvnat6.c

petsc-master 2019-06-15
Report Typos and Errors
  1:  #include <../src/mat/impls/baij/seq/baij.h>
  2:  #include <petsc/private/kernels/blockinvert.h>


  5: PetscErrorCode MatSolve_SeqBAIJ_6_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx)
  6: {
  7:   Mat_SeqBAIJ       *a = (Mat_SeqBAIJ*)A->data;
  8:   PetscInt          i,nz,idx,idt,jdx;
  9:   PetscErrorCode    ierr;
 10:   const PetscInt    *diag = a->diag,*vi,n=a->mbs,*ai=a->i,*aj=a->j;
 11:   const MatScalar   *aa   =a->a,*v;
 12:   PetscScalar       *x,s1,s2,s3,s4,s5,s6,x1,x2,x3,x4,x5,x6;
 13:   const PetscScalar *b;

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

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

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

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

123:     x[idx]   = s1;
124:     x[1+idx] = s2;
125:     x[2+idx] = s3;
126:     x[3+idx] = s4;
127:     x[4+idx] = s5;
128:     x[5+idx] = s6;
129:   }

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

160:   VecRestoreArrayRead(bb,&b);
161:   VecRestoreArray(xx,&x);
162:   PetscLogFlops(2.0*bs2*(a->nz) - bs*A->cmap->n);
163:   return(0);
164: }