#include <../src/mat/impls/baij/seq/baij.h> #include /* Block operations are done by accessing one column at at time */ /* Default MatSolve for block size 11 */ PetscErrorCode MatSolve_SeqBAIJ_11_NaturalOrdering(Mat A,Vec bb,Vec xx) { Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; PetscErrorCode ierr; const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*adiag=a->diag,*vi,bs=A->rmap->bs,bs2=a->bs2; PetscInt i,k,nz,idx,idt,m; const MatScalar *aa=a->a,*v; PetscScalar s[11]; PetscScalar *x,xv; const PetscScalar *b; PetscFunctionBegin; ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); ierr = VecGetArray(xx,&x);CHKERRQ(ierr); /* forward solve the lower triangular */ for (i=0; i=0; i--) { v = aa + bs2*(adiag[i+1]+1); vi = aj + adiag[i+1]+1; nz = adiag[i] - adiag[i+1] - 1; idt = bs*i; s[0] = x[idt]; s[1] = x[1+idt]; s[2] = x[2+idt]; s[3] = x[3+idt]; s[4] = x[4+idt]; 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]; s[10] = x[10+idt]; for (m=0; mnz) - bs*A->cmap->n);CHKERRQ(ierr); PetscFunctionReturn(0); }