1 #include <../src/mat/impls/baij/seq/baij.h> 2 3 PetscErrorCode MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx) 4 { 5 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 6 PetscErrorCode ierr; 7 PetscInt i,nz,idx,idt,oidx; 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 s1,s2,x1,x2,*x; 11 12 PetscFunctionBegin; 13 ierr = VecCopy(bb,xx);CHKERRQ(ierr); 14 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 15 16 /* forward solve the U^T */ 17 idx = 0; 18 for (i=0; i<n; i++) { 19 20 v = aa + 4*diag[i]; 21 /* multiply by the inverse of the block diagonal */ 22 x1 = x[idx]; x2 = x[1+idx]; 23 s1 = v[0]*x1 + v[1]*x2; 24 s2 = v[2]*x1 + v[3]*x2; 25 v += 4; 26 27 vi = aj + diag[i] + 1; 28 nz = ai[i+1] - diag[i] - 1; 29 while (nz--) { 30 oidx = 2*(*vi++); 31 x[oidx] -= v[0]*s1 + v[1]*s2; 32 x[oidx+1] -= v[2]*s1 + v[3]*s2; 33 v += 4; 34 } 35 x[idx] = s1;x[1+idx] = s2; 36 idx += 2; 37 } 38 /* backward solve the L^T */ 39 for (i=n-1; i>=0; i--) { 40 v = aa + 4*diag[i] - 4; 41 vi = aj + diag[i] - 1; 42 nz = diag[i] - ai[i]; 43 idt = 2*i; 44 s1 = x[idt]; s2 = x[1+idt]; 45 while (nz--) { 46 idx = 2*(*vi--); 47 x[idx] -= v[0]*s1 + v[1]*s2; 48 x[idx+1] -= v[2]*s1 + v[3]*s2; 49 v -= 4; 50 } 51 } 52 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 53 ierr = PetscLogFlops(2.0*4.0*(a->nz) - 2.0*A->cmap->n);CHKERRQ(ierr); 54 PetscFunctionReturn(0); 55 } 56 57 PetscErrorCode MatSolveTranspose_SeqBAIJ_2_NaturalOrdering(Mat A,Vec bb,Vec xx) 58 { 59 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 60 PetscErrorCode ierr; 61 const PetscInt n=a->mbs,*vi,*ai=a->i,*aj=a->j,*diag=a->diag; 62 PetscInt nz,idx,idt,j,i,oidx; 63 const PetscInt bs =A->rmap->bs,bs2=a->bs2; 64 const MatScalar *aa=a->a,*v; 65 PetscScalar s1,s2,x1,x2,*x; 66 67 PetscFunctionBegin; 68 ierr = VecCopy(bb,xx);CHKERRQ(ierr); 69 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 70 71 /* forward solve the U^T */ 72 idx = 0; 73 for (i=0; i<n; i++) { 74 v = aa + bs2*diag[i]; 75 /* multiply by the inverse of the block diagonal */ 76 x1 = x[idx]; x2 = x[1+idx]; 77 s1 = v[0]*x1 + v[1]*x2; 78 s2 = v[2]*x1 + v[3]*x2; 79 v -= bs2; 80 81 vi = aj + diag[i] - 1; 82 nz = diag[i] - diag[i+1] - 1; 83 for (j=0; j>-nz; j--) { 84 oidx = bs*vi[j]; 85 x[oidx] -= v[0]*s1 + v[1]*s2; 86 x[oidx+1] -= v[2]*s1 + v[3]*s2; 87 v -= bs2; 88 } 89 x[idx] = s1;x[1+idx] = s2; 90 idx += bs; 91 } 92 /* backward solve the L^T */ 93 for (i=n-1; i>=0; i--) { 94 v = aa + bs2*ai[i]; 95 vi = aj + ai[i]; 96 nz = ai[i+1] - ai[i]; 97 idt = bs*i; 98 s1 = x[idt]; s2 = x[1+idt]; 99 for (j=0; j<nz; j++) { 100 idx = bs*vi[j]; 101 x[idx] -= v[0]*s1 + v[1]*s2; 102 x[idx+1] -= v[2]*s1 + v[3]*s2; 103 v += bs2; 104 } 105 } 106 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 107 ierr = PetscLogFlops(2.0*bs2*(a->nz) - bs*A->cmap->n);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110