2:  #include <../src/mat/utils/freespace.h>

  4: PetscErrorCode PetscFreeSpaceGet(PetscInt n,PetscFreeSpaceList *list)
  5: {
  6:   PetscFreeSpaceList a;
  7:   PetscErrorCode     ierr;

 10:   PetscNew(&a);
 11:   PetscMalloc1(n,&(a->array_head));

 13:   a->array            = a->array_head;
 14:   a->local_remaining  = n;
 15:   a->local_used       = 0;
 16:   a->total_array_size = 0;
 17:   a->more_space       = NULL;

 19:   if (*list) {
 20:     (*list)->more_space = a;
 21:     a->total_array_size = (*list)->total_array_size;
 22:   }

 24:   a->total_array_size += n;
 25:   *list                =  a;
 26:   return(0);
 27: }

 29: PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head,PetscInt *space)
 30: {
 31:   PetscFreeSpaceList a;
 32:   PetscErrorCode     ierr;

 35:   while ((*head)) {
 36:     a      =  (*head)->more_space;
 37:      PetscArraycpy(space,(*head)->array_head,(*head)->local_used);
 38:     space += (*head)->local_used;
 39:      PetscFree((*head)->array_head);
 40:      PetscFree(*head);
 41:     *head  =  a;
 42:   }
 43:   return(0);
 44: }

 46: /*
 47:   PetscFreeSpaceContiguous_LU -
 48:     Copy a linket list obtained from matrix symbolic ILU or LU factorization into a contiguous array
 49:   that enables an efficient matrix triangular solve.

 51:    Input Parameters:
 52: +  head - linked list of column indices obtained from matrix symbolic ILU or LU factorization
 53: .  space - an allocated array with length nnz of factored matrix.
 54: .  n - order of the matrix
 55: .  bi - row pointer of factored matrix L with length n+1.
 56: -  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.

 58:    Output Parameter:
 59: .  space - column indices are copied into this array with contiguous layout of L and U

 61:    See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U
 62: */
 63: PetscErrorCode PetscFreeSpaceContiguous_LU(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *bi,PetscInt *bdiag)
 64: {
 65:   PetscFreeSpaceList a;
 66:   PetscErrorCode     ierr;
 67:   PetscInt           row,nnz,*bj,*array,total,bi_temp;
 68:   PetscInt           nnzL,nnzU;

 71:   bi_temp = bi[n];
 72:   row     = 0;
 73:   total   = 0;
 74:   nnzL    = bdiag[0];
 75:   while ((*head)) {
 76:     total += (*head)->local_used;
 77:     array  = (*head)->array_head;

 79:     while (row < n) {
 80:       if (bi[row+1] > total) break;
 81:       /* copy array entries into bj for this row */
 82:       nnz = bi[row+1] - bi[row];
 83:       /* set bi[row] for new datastruct */
 84:       if (row == 0) {
 85:         bi[row] = 0;
 86:       } else {
 87:         bi[row] = bi[row-1] + nnzL; /* nnzL of previous row */
 88:       }

 90:       /* L part */
 91:       nnzL = bdiag[row];
 92:       bj   = space+bi[row];
 93:       PetscArraycpy(bj,array,nnzL);

 95:       /* diagonal entry */
 96:       bdiag[row]        = bi_temp - 1;
 97:       space[bdiag[row]] = row;

 99:       /* U part */
100:       nnzU    = nnz - nnzL;
101:       bi_temp = bi_temp - nnzU;
102:       nnzU--;       /* exclude diagonal */
103:       bj     = space + bi_temp;
104:       PetscArraycpy(bj,array+nnzL+1,nnzU);
105:       array += nnz;
106:       row++;
107:     }

109:     a     = (*head)->more_space;
110:     PetscFree((*head)->array_head);
111:     PetscFree(*head);
112:     *head = a;
113:   }
114:   if (n) {
115:     bi[n]    = bi[n-1] + nnzL;
116:     bdiag[n] = bdiag[n-1]-1;
117:   }
118:   return(0);
119: }

121: /*
122:   PetscFreeSpaceContiguous_Cholesky -
123:     Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array
124:   that enables an efficient matrix triangular solve.

126:    Input Parameters:
127: +  head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization
128: .  space - an allocated array with length nnz of factored matrix.
129: .  n - order of the matrix
130: .  ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix.
131: -  udiag - array of length n.

133:    Output Parameter:
134: +  space - column indices are copied into this array with contiguous layout of U, with diagonal located as the last entry in each row
135: -  udiag - indices of diagonal entries

137:    See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description.
138: */

140: PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *ui,PetscInt *udiag)
141: {
142:   PetscFreeSpaceList a;
143:   PetscErrorCode     ierr;
144:   PetscInt           row,nnz,*uj,*array,total;

147:   row   = 0;
148:   total = 0;
149:   while (*head) {
150:     total += (*head)->local_used;
151:     array  = (*head)->array_head;

153:     while (row < n) {
154:       if (ui[row+1] > total) break;
155:       udiag[row] = ui[row+1] - 1;     /* points to the last entry of U(row,:) */
156:       nnz        = ui[row+1] - ui[row] - 1; /* exclude diagonal */
157:       uj         = space + ui[row];
158:       PetscArraycpy(uj,array+1,nnz);
159:       uj[nnz]    = array[0]; /* diagonal */
160:       array     += nnz + 1;
161:       row++;
162:     }

164:     a     = (*head)->more_space;
165:     PetscFree((*head)->array_head);
166:     PetscFree(*head);
167:     *head = a;
168:   }
169:   return(0);
170: }

172: PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head)
173: {
174:   PetscFreeSpaceList a;
175:   PetscErrorCode     ierr;

178:   while ((head)) {
179:     a    = (head)->more_space;
180:     PetscFree((head)->array_head);
181:     PetscFree(head);
182:     head = a;
183:   }
184:   return(0);
185: }