1 #include <../src/mat/utils/freespace.h>
2
PetscFreeSpaceGet(PetscInt n,PetscFreeSpaceList * list)3 PetscErrorCode PetscFreeSpaceGet(PetscInt n, PetscFreeSpaceList *list)
4 {
5 PetscFreeSpaceList a;
6
7 PetscFunctionBegin;
8 PetscCall(PetscNew(&a));
9 PetscCall(PetscMalloc1(n, &a->array_head));
10
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;
16
17 if (*list) {
18 (*list)->more_space = a;
19 a->total_array_size = (*list)->total_array_size;
20 }
21
22 a->total_array_size += n;
23 *list = a;
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
PetscFreeSpaceContiguous(PetscFreeSpaceList * head,PetscInt * space)27 PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head, PetscInt *space)
28 {
29 PetscFreeSpaceList a;
30
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 }
42
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.
47
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.
54
55 Output Parameter:
56 . space - column indices are copied into this array with contiguous layout of L and U
57
58 See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U
59 */
PetscFreeSpaceContiguous_LU(PetscFreeSpaceList * head,PetscInt * space,PetscInt n,PetscInt * bi,PetscInt * bdiag)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;
65
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;
74
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 }
85
86 /* L part */
87 nnzL = bdiag[row];
88 bj = space + bi[row];
89 PetscCall(PetscArraycpy(bj, array, nnzL));
90
91 /* diagonal entry */
92 bdiag[row] = bi_temp - 1;
93 space[bdiag[row]] = row;
94
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 }
104
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 }
116
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.
121
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.
128
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
132
133 See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description.
134 */
135
PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList * head,PetscInt * space,PetscInt n,PetscInt * ui,PetscInt * udiag)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;
140
141 PetscFunctionBegin;
142 row = 0;
143 total = 0;
144 while (*head) {
145 total += (*head)->local_used;
146 array = (*head)->array_head;
147
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 }
158
159 a = (*head)->more_space;
160 PetscCall(PetscFree((*head)->array_head));
161 PetscCall(PetscFree(*head));
162 *head = a;
163 }
164 PetscFunctionReturn(PETSC_SUCCESS);
165 }
166
PetscFreeSpaceDestroy(PetscFreeSpaceList head)167 PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head)
168 {
169 PetscFreeSpaceList a;
170
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 }
180