Actual source code: freespace.c
2: #include <../src/mat/utils/freespace.h>
4: PetscErrorCode PetscFreeSpaceGet(PetscInt n, PetscFreeSpaceList *list)
5: {
6: PetscFreeSpaceList a;
8: PetscFunctionBegin;
9: PetscCall(PetscNew(&a));
10: PetscCall(PetscMalloc1(n, &(a->array_head)));
12: a->array = a->array_head;
13: a->local_remaining = n;
14: a->local_used = 0;
15: a->total_array_size = 0;
16: a->more_space = NULL;
18: if (*list) {
19: (*list)->more_space = a;
20: a->total_array_size = (*list)->total_array_size;
21: }
23: a->total_array_size += n;
24: *list = a;
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head, PetscInt *space)
29: {
30: PetscFreeSpaceList a;
32: PetscFunctionBegin;
33: while ((*head)) {
34: a = (*head)->more_space;
35: PetscCall(PetscArraycpy(space, (*head)->array_head, (*head)->local_used));
36: space += (*head)->local_used;
37: PetscCall(PetscFree((*head)->array_head));
38: PetscCall(PetscFree(*head));
39: *head = a;
40: }
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: /*
45: PetscFreeSpaceContiguous_LU -
46: Copy a linket list obtained from matrix symbolic ILU or LU factorization into a contiguous array
47: that enables an efficient matrix triangular solve.
49: Input Parameters:
50: + head - linked list of column indices obtained from matrix symbolic ILU or LU factorization
51: . space - an allocated array with length nnz of factored matrix.
52: . n - order of the matrix
53: . bi - row pointer of factored matrix L with length n+1.
54: - 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.
56: Output Parameter:
57: . space - column indices are copied into this array with contiguous layout of L and U
59: See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U
60: */
61: PetscErrorCode PetscFreeSpaceContiguous_LU(PetscFreeSpaceList *head, PetscInt *space, PetscInt n, PetscInt *bi, PetscInt *bdiag)
62: {
63: PetscFreeSpaceList a;
64: PetscInt row, nnz, *bj, *array, total, bi_temp;
65: PetscInt nnzL, nnzU;
67: PetscFunctionBegin;
68: bi_temp = bi[n];
69: row = 0;
70: total = 0;
71: nnzL = bdiag[0];
72: while ((*head)) {
73: total += (*head)->local_used;
74: array = (*head)->array_head;
76: while (row < n) {
77: if (bi[row + 1] > total) break;
78: /* copy array entries into bj for this row */
79: nnz = bi[row + 1] - bi[row];
80: /* set bi[row] for new datastruct */
81: if (row == 0) {
82: bi[row] = 0;
83: } else {
84: bi[row] = bi[row - 1] + nnzL; /* nnzL of previous row */
85: }
87: /* L part */
88: nnzL = bdiag[row];
89: bj = space + bi[row];
90: PetscCall(PetscArraycpy(bj, array, nnzL));
92: /* diagonal entry */
93: bdiag[row] = bi_temp - 1;
94: space[bdiag[row]] = row;
96: /* U part */
97: nnzU = nnz - nnzL;
98: bi_temp = bi_temp - nnzU;
99: nnzU--; /* exclude diagonal */
100: bj = space + bi_temp;
101: PetscCall(PetscArraycpy(bj, array + nnzL + 1, nnzU));
102: array += nnz;
103: row++;
104: }
106: a = (*head)->more_space;
107: PetscCall(PetscFree((*head)->array_head));
108: PetscCall(PetscFree(*head));
109: *head = a;
110: }
111: if (n) {
112: bi[n] = bi[n - 1] + nnzL;
113: bdiag[n] = bdiag[n - 1] - 1;
114: }
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: /*
119: PetscFreeSpaceContiguous_Cholesky -
120: Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array
121: that enables an efficient matrix triangular solve.
123: Input Parameters:
124: + head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization
125: . space - an allocated array with length nnz of factored matrix.
126: . n - order of the matrix
127: . ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix.
128: - udiag - array of length n.
130: Output Parameter:
131: + space - column indices are copied into this array with contiguous layout of U, with diagonal located as the last entry in each row
132: - udiag - indices of diagonal entries
134: See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description.
135: */
137: PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head, PetscInt *space, PetscInt n, PetscInt *ui, PetscInt *udiag)
138: {
139: PetscFreeSpaceList a;
140: PetscInt row, nnz, *uj, *array, total;
142: PetscFunctionBegin;
143: row = 0;
144: total = 0;
145: while (*head) {
146: total += (*head)->local_used;
147: array = (*head)->array_head;
149: while (row < n) {
150: if (ui[row + 1] > total) break;
151: udiag[row] = ui[row + 1] - 1; /* points to the last entry of U(row,:) */
152: nnz = ui[row + 1] - ui[row] - 1; /* exclude diagonal */
153: uj = space + ui[row];
154: PetscCall(PetscArraycpy(uj, array + 1, nnz));
155: uj[nnz] = array[0]; /* diagonal */
156: array += nnz + 1;
157: row++;
158: }
160: a = (*head)->more_space;
161: PetscCall(PetscFree((*head)->array_head));
162: PetscCall(PetscFree(*head));
163: *head = a;
164: }
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head)
169: {
170: PetscFreeSpaceList a;
172: PetscFunctionBegin;
173: while ((head)) {
174: a = (head)->more_space;
175: PetscCall(PetscFree((head)->array_head));
176: PetscCall(PetscFree(head));
177: head = a;
178: }
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }