1*2c733ed4SBarry Smith #include <../src/mat/impls/baij/seq/baij.h> 2*2c733ed4SBarry Smith 3*2c733ed4SBarry Smith PetscErrorCode MatSolveTranspose_SeqBAIJ_1_NaturalOrdering(Mat A,Vec bb,Vec xx) 4*2c733ed4SBarry Smith { 5*2c733ed4SBarry Smith Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 6*2c733ed4SBarry Smith PetscErrorCode ierr; 7*2c733ed4SBarry Smith const PetscInt *adiag = a->diag,*ai = a->i,*aj = a->j,*vi; 8*2c733ed4SBarry Smith PetscInt i,n = a->mbs,j; 9*2c733ed4SBarry Smith PetscInt nz; 10*2c733ed4SBarry Smith PetscScalar *x,*tmp,s1; 11*2c733ed4SBarry Smith const MatScalar *aa = a->a,*v; 12*2c733ed4SBarry Smith const PetscScalar *b; 13*2c733ed4SBarry Smith 14*2c733ed4SBarry Smith PetscFunctionBegin; 15*2c733ed4SBarry Smith ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 16*2c733ed4SBarry Smith ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 17*2c733ed4SBarry Smith tmp = a->solve_work; 18*2c733ed4SBarry Smith 19*2c733ed4SBarry Smith 20*2c733ed4SBarry Smith /* copy the b into temp work space according to permutation */ 21*2c733ed4SBarry Smith for (i=0; i<n; i++) tmp[i] = b[i]; 22*2c733ed4SBarry Smith 23*2c733ed4SBarry Smith /* forward solve the U^T */ 24*2c733ed4SBarry Smith for (i=0; i<n; i++) { 25*2c733ed4SBarry Smith v = aa + adiag[i+1] + 1; 26*2c733ed4SBarry Smith vi = aj + adiag[i+1] + 1; 27*2c733ed4SBarry Smith nz = adiag[i] - adiag[i+1] - 1; 28*2c733ed4SBarry Smith s1 = tmp[i]; 29*2c733ed4SBarry Smith s1 *= v[nz]; /* multiply by inverse of diagonal entry */ 30*2c733ed4SBarry Smith for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j]; 31*2c733ed4SBarry Smith tmp[i] = s1; 32*2c733ed4SBarry Smith } 33*2c733ed4SBarry Smith 34*2c733ed4SBarry Smith /* backward solve the L^T */ 35*2c733ed4SBarry Smith for (i=n-1; i>=0; i--) { 36*2c733ed4SBarry Smith v = aa + ai[i]; 37*2c733ed4SBarry Smith vi = aj + ai[i]; 38*2c733ed4SBarry Smith nz = ai[i+1] - ai[i]; 39*2c733ed4SBarry Smith s1 = tmp[i]; 40*2c733ed4SBarry Smith for (j=0; j<nz; j++) tmp[vi[j]] -= s1*v[j]; 41*2c733ed4SBarry Smith } 42*2c733ed4SBarry Smith 43*2c733ed4SBarry Smith /* copy tmp into x according to permutation */ 44*2c733ed4SBarry Smith for (i=0; i<n; i++) x[i] = tmp[i]; 45*2c733ed4SBarry Smith 46*2c733ed4SBarry Smith ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 47*2c733ed4SBarry Smith ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 48*2c733ed4SBarry Smith 49*2c733ed4SBarry Smith ierr = PetscLogFlops(2.0*a->nz-A->cmap->n);CHKERRQ(ierr); 50*2c733ed4SBarry Smith PetscFunctionReturn(0); 51*2c733ed4SBarry Smith } 52*2c733ed4SBarry Smith 53*2c733ed4SBarry Smith PetscErrorCode MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx) 54*2c733ed4SBarry Smith { 55*2c733ed4SBarry Smith Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 56*2c733ed4SBarry Smith PetscErrorCode ierr; 57*2c733ed4SBarry Smith PetscInt i,nz; 58*2c733ed4SBarry Smith const PetscInt *diag = a->diag,n=a->mbs,*vi,*ai=a->i,*aj=a->j; 59*2c733ed4SBarry Smith const MatScalar *aa =a->a,*v; 60*2c733ed4SBarry Smith PetscScalar s1,*x; 61*2c733ed4SBarry Smith 62*2c733ed4SBarry Smith PetscFunctionBegin; 63*2c733ed4SBarry Smith ierr = VecCopy(bb,xx);CHKERRQ(ierr); 64*2c733ed4SBarry Smith ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 65*2c733ed4SBarry Smith 66*2c733ed4SBarry Smith /* forward solve the U^T */ 67*2c733ed4SBarry Smith for (i=0; i<n; i++) { 68*2c733ed4SBarry Smith 69*2c733ed4SBarry Smith v = aa + diag[i]; 70*2c733ed4SBarry Smith /* multiply by the inverse of the block diagonal */ 71*2c733ed4SBarry Smith s1 = (*v++)*x[i]; 72*2c733ed4SBarry Smith vi = aj + diag[i] + 1; 73*2c733ed4SBarry Smith nz = ai[i+1] - diag[i] - 1; 74*2c733ed4SBarry Smith while (nz--) { 75*2c733ed4SBarry Smith x[*vi++] -= (*v++)*s1; 76*2c733ed4SBarry Smith } 77*2c733ed4SBarry Smith x[i] = s1; 78*2c733ed4SBarry Smith } 79*2c733ed4SBarry Smith /* backward solve the L^T */ 80*2c733ed4SBarry Smith for (i=n-1; i>=0; i--) { 81*2c733ed4SBarry Smith v = aa + diag[i] - 1; 82*2c733ed4SBarry Smith vi = aj + diag[i] - 1; 83*2c733ed4SBarry Smith nz = diag[i] - ai[i]; 84*2c733ed4SBarry Smith s1 = x[i]; 85*2c733ed4SBarry Smith while (nz--) { 86*2c733ed4SBarry Smith x[*vi--] -= (*v--)*s1; 87*2c733ed4SBarry Smith } 88*2c733ed4SBarry Smith } 89*2c733ed4SBarry Smith ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 90*2c733ed4SBarry Smith ierr = PetscLogFlops(2.0*(a->nz) - A->cmap->n);CHKERRQ(ierr); 91*2c733ed4SBarry Smith PetscFunctionReturn(0); 92*2c733ed4SBarry Smith } 93*2c733ed4SBarry Smith 94