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