Actual source code: baijsolvnat11.c

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

  4: /* Block operations are done by accessing one column at at time */
  5: /* Default MatSolve for block size 11 */

  7: PetscErrorCode MatSolve_SeqBAIJ_11_NaturalOrdering(Mat A,Vec bb,Vec xx)
  8: {
  9:   Mat_SeqBAIJ       *a=(Mat_SeqBAIJ*)A->data;
 10:   const PetscInt    n=a->mbs,*ai=a->i,*aj=a->j,*adiag=a->diag,*vi,bs=A->rmap->bs,bs2=a->bs2;
 11:   PetscInt          i,k,nz,idx,idt,m;
 12:   const MatScalar   *aa=a->a,*v;
 13:   PetscScalar       s[11];
 14:   PetscScalar       *x,xv;
 15:   const PetscScalar *b;

 17:   VecGetArrayRead(bb,&b);
 18:   VecGetArray(xx,&x);

 20:   /* forward solve the lower triangular */
 21:   for (i=0; i<n; i++) {
 22:     v         = aa + bs2*ai[i];
 23:     vi        = aj + ai[i];
 24:     nz        = ai[i+1] - ai[i];
 25:     idt       = bs*i;
 26:     x[idt]    = b[idt];    x[1+idt]  = b[1+idt];  x[2+idt]  = b[2+idt];  x[3+idt]  = b[3+idt];  x[4+idt]  = b[4+idt];
 27:     x[5+idt]  = b[5+idt];  x[6+idt]  = b[6+idt];  x[7+idt]  = b[7+idt];  x[8+idt]  = b[8+idt];  x[9+idt] = b[9+idt];
 28:     x[10+idt] = b[10+idt];
 29:     for (m=0; m<nz; m++) {
 30:       idx = bs*vi[m];
 31:       for (k=0; k<11; k++) {
 32:         xv         = x[k + idx];
 33:         x[idt]    -= v[0]*xv;
 34:         x[1+idt]  -= v[1]*xv;
 35:         x[2+idt]  -= v[2]*xv;
 36:         x[3+idt]  -= v[3]*xv;
 37:         x[4+idt]  -= v[4]*xv;
 38:         x[5+idt]  -= v[5]*xv;
 39:         x[6+idt]  -= v[6]*xv;
 40:         x[7+idt]  -= v[7]*xv;
 41:         x[8+idt]  -= v[8]*xv;
 42:         x[9+idt]  -= v[9]*xv;
 43:         x[10+idt] -= v[10]*xv;
 44:         v         += 11;
 45:       }
 46:     }
 47:   }
 48:   /* backward solve the upper triangular */
 49:   for (i=n-1; i>=0; i--) {
 50:     v     = aa + bs2*(adiag[i+1]+1);
 51:     vi    = aj + adiag[i+1]+1;
 52:     nz    = adiag[i] - adiag[i+1] - 1;
 53:     idt   = bs*i;
 54:     s[0]  = x[idt];    s[1]  = x[1+idt];  s[2]  = x[2+idt];  s[3]  = x[3+idt];  s[4]  = x[4+idt];
 55:     s[5]  = x[5+idt];  s[6]  = x[6+idt];  s[7]  = x[7+idt];  s[8]  = x[8+idt];  s[9]  = x[9+idt];
 56:     s[10] = x[10+idt];

 58:     for (m=0; m<nz; m++) {
 59:       idx = bs*vi[m];
 60:       for (k=0; k<11; k++) {
 61:         xv     = x[k + idx];
 62:         s[0]  -= v[0]*xv;
 63:         s[1]  -= v[1]*xv;
 64:         s[2]  -= v[2]*xv;
 65:         s[3]  -= v[3]*xv;
 66:         s[4]  -= v[4]*xv;
 67:         s[5]  -= v[5]*xv;
 68:         s[6]  -= v[6]*xv;
 69:         s[7]  -= v[7]*xv;
 70:         s[8]  -= v[8]*xv;
 71:         s[9]  -= v[9]*xv;
 72:         s[10] -= v[10]*xv;
 73:         v     += 11;
 74:       }
 75:     }
 76:     PetscArrayzero(x+idt,bs);
 77:     for (k=0; k<11; k++) {
 78:       x[idt]    += v[0]*s[k];
 79:       x[1+idt]  += v[1]*s[k];
 80:       x[2+idt]  += v[2]*s[k];
 81:       x[3+idt]  += v[3]*s[k];
 82:       x[4+idt]  += v[4]*s[k];
 83:       x[5+idt]  += v[5]*s[k];
 84:       x[6+idt]  += v[6]*s[k];
 85:       x[7+idt]  += v[7]*s[k];
 86:       x[8+idt]  += v[8]*s[k];
 87:       x[9+idt]  += v[9]*s[k];
 88:       x[10+idt] += v[10]*s[k];
 89:       v         += 11;
 90:     }
 91:   }
 92:   VecRestoreArrayRead(bb,&b);
 93:   VecRestoreArray(xx,&x);
 94:   PetscLogFlops(2.0*bs2*(a->nz) - bs*A->cmap->n);
 95:   return 0;
 96: }