12c733ed4SBarry Smith #include <../src/mat/impls/baij/seq/baij.h> 22c733ed4SBarry Smith 39371c9d4SSatish Balay PetscErrorCode MatSolveTranspose_SeqBAIJ_1_NaturalOrdering(Mat A, Vec bb, Vec xx) { 42c733ed4SBarry Smith Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 52c733ed4SBarry Smith const PetscInt *adiag = a->diag, *ai = a->i, *aj = a->j, *vi; 62c733ed4SBarry Smith PetscInt i, n = a->mbs, j; 72c733ed4SBarry Smith PetscInt nz; 82c733ed4SBarry Smith PetscScalar *x, *tmp, s1; 92c733ed4SBarry Smith const MatScalar *aa = a->a, *v; 102c733ed4SBarry Smith const PetscScalar *b; 112c733ed4SBarry Smith 122c733ed4SBarry Smith PetscFunctionBegin; 139566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(bb, &b)); 149566063dSJacob Faibussowitsch PetscCall(VecGetArray(xx, &x)); 152c733ed4SBarry Smith tmp = a->solve_work; 162c733ed4SBarry Smith 172c733ed4SBarry Smith /* copy the b into temp work space according to permutation */ 182c733ed4SBarry Smith for (i = 0; i < n; i++) tmp[i] = b[i]; 192c733ed4SBarry Smith 202c733ed4SBarry Smith /* forward solve the U^T */ 212c733ed4SBarry Smith for (i = 0; i < n; i++) { 222c733ed4SBarry Smith v = aa + adiag[i + 1] + 1; 232c733ed4SBarry Smith vi = aj + adiag[i + 1] + 1; 242c733ed4SBarry Smith nz = adiag[i] - adiag[i + 1] - 1; 252c733ed4SBarry Smith s1 = tmp[i]; 262c733ed4SBarry Smith s1 *= v[nz]; /* multiply by inverse of diagonal entry */ 272c733ed4SBarry Smith for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j]; 282c733ed4SBarry Smith tmp[i] = s1; 292c733ed4SBarry Smith } 302c733ed4SBarry Smith 312c733ed4SBarry Smith /* backward solve the L^T */ 322c733ed4SBarry Smith for (i = n - 1; i >= 0; i--) { 332c733ed4SBarry Smith v = aa + ai[i]; 342c733ed4SBarry Smith vi = aj + ai[i]; 352c733ed4SBarry Smith nz = ai[i + 1] - ai[i]; 362c733ed4SBarry Smith s1 = tmp[i]; 372c733ed4SBarry Smith for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j]; 382c733ed4SBarry Smith } 392c733ed4SBarry Smith 402c733ed4SBarry Smith /* copy tmp into x according to permutation */ 412c733ed4SBarry Smith for (i = 0; i < n; i++) x[i] = tmp[i]; 422c733ed4SBarry Smith 439566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(bb, &b)); 449566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(xx, &x)); 452c733ed4SBarry Smith 469566063dSJacob Faibussowitsch PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n)); 472c733ed4SBarry Smith PetscFunctionReturn(0); 482c733ed4SBarry Smith } 492c733ed4SBarry Smith 509371c9d4SSatish Balay PetscErrorCode MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace(Mat A, Vec bb, Vec xx) { 512c733ed4SBarry Smith Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 522c733ed4SBarry Smith PetscInt i, nz; 532c733ed4SBarry Smith const PetscInt *diag = a->diag, n = a->mbs, *vi, *ai = a->i, *aj = a->j; 542c733ed4SBarry Smith const MatScalar *aa = a->a, *v; 552c733ed4SBarry Smith PetscScalar s1, *x; 562c733ed4SBarry Smith 572c733ed4SBarry Smith PetscFunctionBegin; 589566063dSJacob Faibussowitsch PetscCall(VecCopy(bb, xx)); 599566063dSJacob Faibussowitsch PetscCall(VecGetArray(xx, &x)); 602c733ed4SBarry Smith 612c733ed4SBarry Smith /* forward solve the U^T */ 622c733ed4SBarry Smith for (i = 0; i < n; i++) { 632c733ed4SBarry Smith v = aa + diag[i]; 642c733ed4SBarry Smith /* multiply by the inverse of the block diagonal */ 652c733ed4SBarry Smith s1 = (*v++) * x[i]; 662c733ed4SBarry Smith vi = aj + diag[i] + 1; 672c733ed4SBarry Smith nz = ai[i + 1] - diag[i] - 1; 68*ad540459SPierre Jolivet while (nz--) x[*vi++] -= (*v++) * s1; 692c733ed4SBarry Smith x[i] = s1; 702c733ed4SBarry Smith } 712c733ed4SBarry Smith /* backward solve the L^T */ 722c733ed4SBarry Smith for (i = n - 1; i >= 0; i--) { 732c733ed4SBarry Smith v = aa + diag[i] - 1; 742c733ed4SBarry Smith vi = aj + diag[i] - 1; 752c733ed4SBarry Smith nz = diag[i] - ai[i]; 762c733ed4SBarry Smith s1 = x[i]; 77*ad540459SPierre Jolivet while (nz--) x[*vi--] -= (*v--) * s1; 782c733ed4SBarry Smith } 799566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(xx, &x)); 809566063dSJacob Faibussowitsch PetscCall(PetscLogFlops(2.0 * (a->nz) - A->cmap->n)); 812c733ed4SBarry Smith PetscFunctionReturn(0); 822c733ed4SBarry Smith } 83