1: #include <../src/mat/utils/freespace.h>

  3: PetscErrorCode PetscFreeSpaceGet(PetscInt n, PetscFreeSpaceList *list)
  4: {
  5:   PetscFreeSpaceList a;

  7:   PetscFunctionBegin;
  8:   PetscCall(PetscNew(&a));
  9:   PetscCall(PetscMalloc1(n, &a->array_head));

 11:   a->array            = a->array_head;
 12:   a->local_remaining  = n;
 13:   a->local_used       = 0;
 14:   a->total_array_size = 0;
 15:   a->more_space       = NULL;

 17:   if (*list) {
 18:     (*list)->more_space = a;
 19:     a->total_array_size = (*list)->total_array_size;
 20:   }

 22:   a->total_array_size += n;
 23:   *list = a;
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

 27: PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head, PetscInt *space)
 28: {
 29:   PetscFreeSpaceList a;

 31:   PetscFunctionBegin;
 32:   while (*head) {
 33:     a = (*head)->more_space;
 34:     PetscCall(PetscArraycpy(space, (*head)->array_head, (*head)->local_used));
 35:     space = PetscSafePointerPlusOffset(space, (*head)->local_used);
 36:     PetscCall(PetscFree((*head)->array_head));
 37:     PetscCall(PetscFree(*head));
 38:     *head = a;
 39:   }
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: /*
 44:   PetscFreeSpaceContiguous_LU -
 45:     Copy a linket list obtained from matrix symbolic ILU or LU factorization into a contiguous array
 46:   that enables an efficient matrix triangular solve.

 48:    Input Parameters:
 49: +  head - linked list of column indices obtained from matrix symbolic ILU or LU factorization
 50: .  space - an allocated array with length nnz of factored matrix.
 51: .  n - order of the matrix
 52: .  bi - row pointer of factored matrix L with length n+1.
 53: -  bdiag - array of length n+1. bdiag[i] points to diagonal of U(i,:), and bdiag[n] points to entry of U(n-1,0)-1.

 55:    Output Parameter:
 56: .  space - column indices are copied into this array with contiguous layout of L and U

 58:    See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U
 59: */
 60: PetscErrorCode PetscFreeSpaceContiguous_LU(PetscFreeSpaceList *head, PetscInt *space, PetscInt n, PetscInt *bi, PetscInt *bdiag)
 61: {
 62:   PetscFreeSpaceList a;
 63:   PetscInt           row, nnz, *bj, *array, total, bi_temp;
 64:   PetscInt           nnzL, nnzU;

 66:   PetscFunctionBegin;
 67:   bi_temp = bi[n];
 68:   row     = 0;
 69:   total   = 0;
 70:   nnzL    = bdiag[0];
 71:   while (*head) {
 72:     total += (*head)->local_used;
 73:     array = (*head)->array_head;

 75:     while (row < n) {
 76:       if (bi[row + 1] > total) break;
 77:       /* copy array entries into bj for this row */
 78:       nnz = bi[row + 1] - bi[row];
 79:       /* set bi[row] for new datastruct */
 80:       if (row == 0) {
 81:         bi[row] = 0;
 82:       } else {
 83:         bi[row] = bi[row - 1] + nnzL; /* nnzL of previous row */
 84:       }

 86:       /* L part */
 87:       nnzL = bdiag[row];
 88:       bj   = space + bi[row];
 89:       PetscCall(PetscArraycpy(bj, array, nnzL));

 91:       /* diagonal entry */
 92:       bdiag[row]        = bi_temp - 1;
 93:       space[bdiag[row]] = row;

 95:       /* U part */
 96:       nnzU    = nnz - nnzL;
 97:       bi_temp = bi_temp - nnzU;
 98:       nnzU--; /* exclude diagonal */
 99:       bj = space + bi_temp;
100:       PetscCall(PetscArraycpy(bj, array + nnzL + 1, nnzU));
101:       array += nnz;
102:       row++;
103:     }

105:     a = (*head)->more_space;
106:     PetscCall(PetscFree((*head)->array_head));
107:     PetscCall(PetscFree(*head));
108:     *head = a;
109:   }
110:   if (n) {
111:     bi[n]    = bi[n - 1] + nnzL;
112:     bdiag[n] = bdiag[n - 1] - 1;
113:   }
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: /*
118:   PetscFreeSpaceContiguous_Cholesky -
119:     Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array
120:   that enables an efficient matrix triangular solve.

122:    Input Parameters:
123: +  head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization
124: .  space - an allocated array with length nnz of factored matrix.
125: .  n - order of the matrix
126: .  ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix.
127: -  udiag - array of length n.

129:    Output Parameter:
130: +  space - column indices are copied into this array with contiguous layout of U, with diagonal located as the last entry in each row
131: -  udiag - indices of diagonal entries

133:    See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description.
134: */

136: PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head, PetscInt *space, PetscInt n, PetscInt *ui, PetscInt *udiag)
137: {
138:   PetscFreeSpaceList a;
139:   PetscInt           row, nnz, *uj, *array, total;

141:   PetscFunctionBegin;
142:   row   = 0;
143:   total = 0;
144:   while (*head) {
145:     total += (*head)->local_used;
146:     array = (*head)->array_head;

148:     while (row < n) {
149:       if (ui[row + 1] > total) break;
150:       udiag[row] = ui[row + 1] - 1;           /* points to the last entry of U(row,:) */
151:       nnz        = ui[row + 1] - ui[row] - 1; /* exclude diagonal */
152:       uj         = space + ui[row];
153:       PetscCall(PetscArraycpy(uj, array + 1, nnz));
154:       uj[nnz] = array[0]; /* diagonal */
155:       array += nnz + 1;
156:       row++;
157:     }

159:     a = (*head)->more_space;
160:     PetscCall(PetscFree((*head)->array_head));
161:     PetscCall(PetscFree(*head));
162:     *head = a;
163:   }
164:   PetscFunctionReturn(PETSC_SUCCESS);
165: }

167: PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head)
168: {
169:   PetscFreeSpaceList a;

171:   PetscFunctionBegin;
172:   while (head) {
173:     a = (head)->more_space;
174:     PetscCall(PetscFree((head)->array_head));
175:     PetscCall(PetscFree(head));
176:     head = a;
177:   }
178:   PetscFunctionReturn(PETSC_SUCCESS);
179: }