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: }