Actual source code: baijsolvtrannat2.c

  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:   VecCopy(bb, xx);
 12:   VecGetArray(xx, &x);

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

 25:     vi = aj + diag[i] + 1;
 26:     nz = ai[i + 1] - diag[i] - 1;
 27:     while (nz--) {
 28:       oidx = 2 * (*vi++);
 29:       x[oidx] -= v[0] * s1 + v[1] * s2;
 30:       x[oidx + 1] -= v[2] * s1 + v[3] * s2;
 31:       v += 4;
 32:     }
 33:     x[idx]     = s1;
 34:     x[1 + idx] = s2;
 35:     idx += 2;
 36:   }
 37:   /* backward solve the L^T */
 38:   for (i = n - 1; i >= 0; i--) {
 39:     v   = aa + 4 * diag[i] - 4;
 40:     vi  = aj + diag[i] - 1;
 41:     nz  = diag[i] - ai[i];
 42:     idt = 2 * i;
 43:     s1  = x[idt];
 44:     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:   VecRestoreArray(xx, &x);
 53:   PetscLogFlops(2.0 * 4.0 * (a->nz) - 2.0 * A->cmap->n);
 54:   return 0;
 55: }

 57: PetscErrorCode MatSolveTranspose_SeqBAIJ_2_NaturalOrdering(Mat A, Vec bb, Vec xx)
 58: {
 59:   Mat_SeqBAIJ     *a = (Mat_SeqBAIJ *)A->data;
 60:   const PetscInt   n = a->mbs, *vi, *ai = a->i, *aj = a->j, *diag = a->diag;
 61:   PetscInt         nz, idx, idt, j, i, oidx;
 62:   const PetscInt   bs = A->rmap->bs, bs2 = a->bs2;
 63:   const MatScalar *aa = a->a, *v;
 64:   PetscScalar      s1, s2, x1, x2, *x;

 66:   VecCopy(bb, xx);
 67:   VecGetArray(xx, &x);

 69:   /* forward solve the U^T */
 70:   idx = 0;
 71:   for (i = 0; i < n; i++) {
 72:     v = aa + bs2 * diag[i];
 73:     /* multiply by the inverse of the block diagonal */
 74:     x1 = x[idx];
 75:     x2 = x[1 + idx];
 76:     s1 = v[0] * x1 + v[1] * x2;
 77:     s2 = v[2] * x1 + v[3] * x2;
 78:     v -= bs2;

 80:     vi = aj + diag[i] - 1;
 81:     nz = diag[i] - diag[i + 1] - 1;
 82:     for (j = 0; j > -nz; j--) {
 83:       oidx = bs * vi[j];
 84:       x[oidx] -= v[0] * s1 + v[1] * s2;
 85:       x[oidx + 1] -= v[2] * s1 + v[3] * s2;
 86:       v -= bs2;
 87:     }
 88:     x[idx]     = s1;
 89:     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];
 99:     s2  = x[1 + idt];
100:     for (j = 0; j < nz; j++) {
101:       idx = bs * vi[j];
102:       x[idx] -= v[0] * s1 + v[1] * s2;
103:       x[idx + 1] -= v[2] * s1 + v[3] * s2;
104:       v += bs2;
105:     }
106:   }
107:   VecRestoreArray(xx, &x);
108:   PetscLogFlops(2.0 * bs2 * (a->nz) - bs * A->cmap->n);
109:   return 0;
110: }