Actual source code: freespace.c
2: #include <../src/mat/utils/freespace.h>
4: PetscErrorCode PetscFreeSpaceGet(PetscInt n, PetscFreeSpaceList *list)
5: {
6: PetscFreeSpaceList a;
8: PetscNew(&a);
9: 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: return 0;
25: }
27: PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head, PetscInt *space)
28: {
29: PetscFreeSpaceList a;
31: while ((*head)) {
32: a = (*head)->more_space;
33: PetscArraycpy(space, (*head)->array_head, (*head)->local_used);
34: space += (*head)->local_used;
35: PetscFree((*head)->array_head);
36: PetscFree(*head);
37: *head = a;
38: }
39: return 0;
40: }
42: /*
43: PetscFreeSpaceContiguous_LU -
44: Copy a linket list obtained from matrix symbolic ILU or LU factorization into a contiguous array
45: that enables an efficient matrix triangular solve.
47: Input Parameters:
48: + head - linked list of column indices obtained from matrix symbolic ILU or LU factorization
49: . space - an allocated array with length nnz of factored matrix.
50: . n - order of the matrix
51: . bi - row pointer of factored matrix L with length n+1.
52: - 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: Output Parameter:
55: . space - column indices are copied into this array with contiguous layout of L and U
57: See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U
58: */
59: PetscErrorCode PetscFreeSpaceContiguous_LU(PetscFreeSpaceList *head, PetscInt *space, PetscInt n, PetscInt *bi, PetscInt *bdiag)
60: {
61: PetscFreeSpaceList a;
62: PetscInt row, nnz, *bj, *array, total, bi_temp;
63: PetscInt nnzL, nnzU;
65: bi_temp = bi[n];
66: row = 0;
67: total = 0;
68: nnzL = bdiag[0];
69: while ((*head)) {
70: total += (*head)->local_used;
71: array = (*head)->array_head;
73: while (row < n) {
74: if (bi[row + 1] > total) break;
75: /* copy array entries into bj for this row */
76: nnz = bi[row + 1] - bi[row];
77: /* set bi[row] for new datastruct */
78: if (row == 0) {
79: bi[row] = 0;
80: } else {
81: bi[row] = bi[row - 1] + nnzL; /* nnzL of previous row */
82: }
84: /* L part */
85: nnzL = bdiag[row];
86: bj = space + bi[row];
87: PetscArraycpy(bj, array, nnzL);
89: /* diagonal entry */
90: bdiag[row] = bi_temp - 1;
91: space[bdiag[row]] = row;
93: /* U part */
94: nnzU = nnz - nnzL;
95: bi_temp = bi_temp - nnzU;
96: nnzU--; /* exclude diagonal */
97: bj = space + bi_temp;
98: PetscArraycpy(bj, array + nnzL + 1, nnzU);
99: array += nnz;
100: row++;
101: }
103: a = (*head)->more_space;
104: PetscFree((*head)->array_head);
105: PetscFree(*head);
106: *head = a;
107: }
108: if (n) {
109: bi[n] = bi[n - 1] + nnzL;
110: bdiag[n] = bdiag[n - 1] - 1;
111: }
112: return 0;
113: }
115: /*
116: PetscFreeSpaceContiguous_Cholesky -
117: Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array
118: that enables an efficient matrix triangular solve.
120: Input Parameters:
121: + head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization
122: . space - an allocated array with length nnz of factored matrix.
123: . n - order of the matrix
124: . ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix.
125: - udiag - array of length n.
127: Output Parameter:
128: + space - column indices are copied into this array with contiguous layout of U, with diagonal located as the last entry in each row
129: - udiag - indices of diagonal entries
131: See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description.
132: */
134: PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head, PetscInt *space, PetscInt n, PetscInt *ui, PetscInt *udiag)
135: {
136: PetscFreeSpaceList a;
137: PetscInt row, nnz, *uj, *array, total;
139: row = 0;
140: total = 0;
141: while (*head) {
142: total += (*head)->local_used;
143: array = (*head)->array_head;
145: while (row < n) {
146: if (ui[row + 1] > total) break;
147: udiag[row] = ui[row + 1] - 1; /* points to the last entry of U(row,:) */
148: nnz = ui[row + 1] - ui[row] - 1; /* exclude diagonal */
149: uj = space + ui[row];
150: PetscArraycpy(uj, array + 1, nnz);
151: uj[nnz] = array[0]; /* diagonal */
152: array += nnz + 1;
153: row++;
154: }
156: a = (*head)->more_space;
157: PetscFree((*head)->array_head);
158: PetscFree(*head);
159: *head = a;
160: }
161: return 0;
162: }
164: PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head)
165: {
166: PetscFreeSpaceList a;
168: while ((head)) {
169: a = (head)->more_space;
170: PetscFree((head)->array_head);
171: PetscFree(head);
172: head = a;
173: }
174: return 0;
175: }