1: #include <../src/mat/impls/baij/seq/baij.h>

  3: PetscErrorCode MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace(Mat A, Vec bb, Vec xx)
  4: {
  5:   Mat_SeqBAIJ     *a = (Mat_SeqBAIJ *)A->data;
  6:   PetscInt         i, nz, idx, idt, oidx;
  7:   const PetscInt  *diag = a->diag, *vi, n = a->mbs, *ai = a->i, *aj = a->j;
  8:   const MatScalar *aa = a->a, *v;
  9:   PetscScalar      s1, s2, x1, x2, *x;

 11:   PetscFunctionBegin;
 12:   PetscCall(VecCopy(bb, xx));
 13:   PetscCall(VecGetArray(xx, &x));

 15:   /* forward solve the U^T */
 16:   idx = 0;
 17:   for (i = 0; i < n; i++) {
 18:     v = aa + 4 * diag[i];
 19:     /* multiply by the inverse of the block diagonal */
 20:     x1 = x[idx];
 21:     x2 = x[1 + idx];
 22:     s1 = v[0] * x1 + v[1] * x2;
 23:     s2 = v[2] * x1 + v[3] * x2;
 24:     v += 4;

 26:     vi = aj + diag[i] + 1;
 27:     nz = ai[i + 1] - diag[i] - 1;
 28:     while (nz--) {
 29:       oidx = 2 * (*vi++);
 30:       x[oidx] -= v[0] * s1 + v[1] * s2;
 31:       x[oidx + 1] -= v[2] * s1 + v[3] * s2;
 32:       v += 4;
 33:     }
 34:     x[idx]     = s1;
 35:     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];
 45:     s2  = x[1 + idt];
 46:     while (nz--) {
 47:       idx = 2 * (*vi--);
 48:       x[idx] -= v[0] * s1 + v[1] * s2;
 49:       x[idx + 1] -= v[2] * s1 + v[3] * s2;
 50:       v -= 4;
 51:     }
 52:   }
 53:   PetscCall(VecRestoreArray(xx, &x));
 54:   PetscCall(PetscLogFlops(2.0 * 4.0 * (a->nz) - 2.0 * A->cmap->n));
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

 58: PetscErrorCode MatSolveTranspose_SeqBAIJ_2_NaturalOrdering(Mat A, Vec bb, Vec xx)
 59: {
 60:   Mat_SeqBAIJ     *a = (Mat_SeqBAIJ *)A->data;
 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;

 67:   PetscFunctionBegin;
 68:   PetscCall(VecCopy(bb, xx));
 69:   PetscCall(VecGetArray(xx, &x));

 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];
 77:     x2 = x[1 + idx];
 78:     s1 = v[0] * x1 + v[1] * x2;
 79:     s2 = v[2] * x1 + v[3] * x2;
 80:     v -= bs2;

 82:     vi = aj + diag[i] - 1;
 83:     nz = diag[i] - diag[i + 1] - 1;
 84:     for (j = 0; j > -nz; j--) {
 85:       oidx = bs * vi[j];
 86:       x[oidx] -= v[0] * s1 + v[1] * s2;
 87:       x[oidx + 1] -= v[2] * s1 + v[3] * s2;
 88:       v -= bs2;
 89:     }
 90:     x[idx]     = s1;
 91:     x[1 + idx] = s2;
 92:     idx += bs;
 93:   }
 94:   /* backward solve the L^T */
 95:   for (i = n - 1; i >= 0; i--) {
 96:     v   = aa + bs2 * ai[i];
 97:     vi  = aj + ai[i];
 98:     nz  = ai[i + 1] - ai[i];
 99:     idt = bs * i;
100:     s1  = x[idt];
101:     s2  = x[1 + idt];
102:     for (j = 0; j < nz; j++) {
103:       idx = bs * vi[j];
104:       x[idx] -= v[0] * s1 + v[1] * s2;
105:       x[idx + 1] -= v[2] * s1 + v[3] * s2;
106:       v += bs2;
107:     }
108:   }
109:   PetscCall(VecRestoreArray(xx, &x));
110:   PetscCall(PetscLogFlops(2.0 * bs2 * (a->nz) - bs * A->cmap->n));
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }